# Symmetry classes of classical stochastic processes

Lucas Sá, Pedro Ribeiro, Tomaž Prosen, Denis Bernard

► To cite this version:

Lucas Sá, Pedro Ribeiro, Tomaž Prosen, Denis Bernard. Symmetry classes of classical stochastic processes. Journal of Statistical Physics, 2025, 192 (3), pp.41. 10.1007/s10955-025-03423-y. hal-04650458

HAL Id: hal-04650458

https://hal.science/hal-04650458v1

Submitted on 26 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Symmetry Classes of Classical Stochastic Processes

Lucas Sá¹ ID · Pedro Ribeiro²,³ ID · Tomaž Prosen⁴,⁵ ID · Denis Bernard⁶ ID

Received: 16 July 2024 / Accepted: 19 February 2025 / Published online: 12 March 2025
© The Author(s) 2025

## Abstract

We perform a systematic symmetry classification of the Markov generators of classical stochastic processes. Our classification scheme is based on the action of involutive symmetry transformations of a real Markov generator, extending the Bernard-LeClair scheme to the arena of classical stochastic processes and leading to a set of up to fifteen allowed symmetry classes. We construct families of solutions of arbitrary matrix dimensions for five of these classes with a simple physical interpretation of particles hopping on multipartite graphs. In the remaining classes, such a simple construction is prevented by the positivity of entries of the generator particular to classical stochastic processes, which imposes a further requirement beyond the usual symmetry classification constraints. We partially overcome this difficulty by resorting to a stochastic optimization algorithm, finding specific examples of generators of small matrix dimensions in six further classes, leaving the existence of the final four allowed classes an open problem. Our symmetry-based results unveil new possibilities in the dynamics of classical stochastic processes: Kramers degeneracy of eigenvalue pairs, dihedral symmetry of the spectra of Markov generators, and time reversal properties of stochastic trajectories and correlation functions.

**Keywords** Stochastic processes · Symmetry classes · Markov generator · Non-Hermitian physics

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Fifteenfold Symmetry Classification</td><td>3</td></tr><tr><td>3</td><td>Analytical Families of Markov Matrices</td><td>6</td></tr><tr><td>3.1</td><td>Class AI, S<sup>2</sup> = R<sub>+</sub>R<sub>+</sub><sup>-T</sup> = R<sub>-</sub>R<sub>-</sub><sup>-T</sup> = 0</td><td>8</td></tr><tr><td>3.2</td><td>Class AI+, S<sup>2</sup> = +1, R<sub>+</sub>R<sub>+</sub><sup>-T</sup> = R<sub>-</sub>R<sub>-</sub><sup>-T</sup> = 0</td><td>8</td></tr><tr><td>3.3</td><td>Class BDI<sup>†</sup>, R<sub>+</sub>R<sub>+</sub><sup>-T</sup> = +1, R<sub>-</sub>R<sub>-</sub><sup>-T</sup> = S<sup>2</sup> = 0</td><td>8</td></tr><tr><td>3.4</td><td>Class CI<sub>+-</sub>, R<sub>+</sub>R<sub>+</sub><sup>-T</sup> = S<sup>2</sup> = +1, R<sub>-</sub>R<sub>-</sub><sup>-T</sup> = -1</td><td>9</td></tr><tr><td>3.5</td><td>Class BDI<sub>++</sub>, R<sub>+</sub>R<sub>+</sub><sup>-T</sup> = R<sub>-</sub>R<sub>-</sub><sup>-T</sup> = S<sup>2</sup> = +1</td><td>9</td></tr><tr><td>3.6</td><td>Absence of Analytic Examples for the Other Classes</td><td>10</td></tr><tr><td>4</td><td>Systematic Numerical Solutions of the Symmetry Classes</td><td>10</td></tr><tr><td>4.1</td><td>Parametrization of the Symmetry Operators</td><td>10</td></tr><tr><td>4.2</td><td>Solution of the Symmetry Constraints</td><td>12</td></tr><tr><td>4.3</td><td>Monte Carlo Sampling</td><td>12</td></tr></table>

Communicated by Christian Maes.

Extended author information available on the last page of the article
<table><tbody><tr><td>4.4</td><td>Results</td><td>13</td></tr><tr><td>5</td><td>Examples and Physical Consequences</td><td>14</td></tr><tr><td>5.1</td><td>Example: Bipartite Graph</td><td>14</td></tr><tr><td>5.2</td><td>Dihedral Symmetry of the Spectrum</td><td>15</td></tr><tr><td>5.3</td><td>Kramers Degeneracy</td><td>16</td></tr><tr><td>5.4</td><td>Detailed Balance and R<sub>+</sub> Symmetry</td><td>17</td></tr><tr><td>5.5</td><td>Time Reversal of Stochastic Trajectories</td><td>17</td></tr><tr><td>5.6</td><td>Time Reversal of Correlation Functions</td><td>18</td></tr><tr><td>6</td><td>Conclusions and Outlook</td><td>18</td></tr><tr><td></td><td>Appendix A Steady-State Blocks of Reducible Markov Matrices</td><td>19</td></tr><tr><td></td><td>A.1 Example: SU(2) Symmetry</td><td>20</td></tr><tr><td></td><td>Appendix B Explicit Examples of Markov Matrices in Different Classes</td><td>22</td></tr><tr><td></td><td>B.1 Class AI</td><td>22</td></tr><tr><td></td><td>B.2 Class DIII<sup>†</sup></td><td>22</td></tr><tr><td></td><td>B.3 Class CI</td><td>23</td></tr><tr><td></td><td>B.4 Class CI++</td><td>23</td></tr><tr><td></td><td>B.5 Class CI</td><td>23</td></tr><tr><td></td><td>B.6 Class CI<sub>-</sub>+</td><td>24</td></tr><tr><td></td><td>References</td><td>24</td></tr></tbody></table>

# 1 Introduction

Random matrix theory pervades modern physics and mathematics [1]. One of its most appealing features is that spectral correlations of complex systems depend only on a small number of symmetries and coincide with those of a random matrix with the same symmetries. Consequently, organizing systems into classes labeled by their symmetries provides extensive information from minimal input. Famous examples include the degree of repulsion between energy levels, the double degeneracy of half-integer-spin time-reversal-symmetric systems, or the universality of conductance fluctuations. These are often the only analytical statements one can make for systems that are not solvable. Quantum mechanics still provides not only most of the applications of random matrix theory but also the most common framework to classify random matrices. Indeed, work on chiral fermions in QCD [2, 3] and disordered superconductors [4] established a tenfold classification of Hermitian Hamiltonians (or, more generally, random matrices) based on unitary and antiunitary symmetries (time-reversal, particle-hole, and chiral), the Altland-Zirnbauer classes. More recently, symmetry classifications of non-Hermitian quantum dynamical generators have also been put forward, for the cases of general non-Hermitian Hamiltonians [5–7], PT-symmetric Hamiltonians [8], and Lindbladians [9–13].

Despite the power and generality of the symmetry classification approach, it has remained mostly circumscribed to the quantum realm. A natural question is, then, how to extend it to classical settings—in particular, nonequilibrium ones described by non-Hermitian matrices— raising the possibility of previously unidentified classes of classical stochastic dynamics, with unique spectral and dynamical consequences. We will focus our attention on classical Markov processes, which describe systems exhibiting memoryless and probabilistic transitions between states, with applications ranging from physical, chemical, and biological systems, to economics. The solution of a classical Markov process is obtained by specifying the time-evolution of a (real) probability vector $p(\sum_{i=1}^{N} p_i = 1)$ of an $N$-dimensional system, governed by the Markov equation

$$ \partial_t p = L p. \qquad (1) $$
Conservation of probability requires that the Markov generator, an $N \times N$ matrix $L$, satisfies the three properties:

$$L_{ij} \in \mathbb{R}, \qquad (2)$$

$$L_{ij} \ge 0, \quad i \neq j, \qquad (3)$$

$$\sum_{i=1}^{N} L_{ij} = 0, \quad \forall j. \qquad (4)$$

Such processes are a natural candidate to investigate symmetry classifications in classical systems, as they can be obtained as a limit (diagonal projection) of quantum Lindbladian evolution, for which the symmetry classification is well-established [9].

In this paper, we show that indeed symmetry classifications are not restricted to the quantum realm, and can be extended to classical stochastic processes. In particular, we show that considering (anti)unitary symmetries is not crucial to developing the classification, thus extending the Altland-Zirnbauer [4] and Bernard-LeClair schemes [5, 6] to general involutive symmetries. We then apply our classificatory scheme to the Markov generators of classical stochastic processes, finding that they can belong at most to fifteen different symmetry classes (Sect. 2). Interestingly, a certain constraint of the quantum symmetry classification that prevents classes with Kramers degeneracy (the exact double degeneracy of every eigenvalue) [9] is absent in the classical limit, illustrating a rich interplay of symmetry and classical stochastic evolution going beyond a mere projecting of quantum Markov dynamics to the classical subspace. Notwithstanding, the classical setting presents its own set of difficulties for the classification. Indeed, contrary to the quantum case, the existence of a preferred basis (where probabilities are real and positive) is at odds with the basis-independent statements sought after by symmetry classifications. Specifically, while we can explicitly build parametric families of Markov generators belonging to five of the fifteen classes, the aforementioned tension prevents the construction of examples in the remaining ten classes if working in the computational basis (Sect. 3). To overcome this difficulty, we developed a stochastic optimization algorithm, which allows us to find examples in six further classes (Sect. 4). The status of the other four classes remains inconclusive, as we could neither build an example for any of them nor exclude its existence from our general fifteenfold classification. Finally, we also discuss consequences for the spectrum of the Markov matrix and for the dynamics, such as Kramers degeneracy of pairs of eigenvalues, dihedral symmetry of the spectrum, time reversal of trajectories and correlation functions, and the role of detailed balance (Sect. 5).

## 2 Fifteenfold Symmetry Classification

In quantum mechanics, to perform the symmetry classification of a non-Hermitian matrix, we look for the antiunitary and unitary symmetry involutions of its irreducible blocks. In a classical stochastic process, there is, a priori, no reason to constrain the symmetries to be unitary or antiunitary. Moreover, the Markov generator $L$ is automatically real and, without any reason to introduce complex numbers at this stage, we therefore demand that the operators implementing the symmetries are themselves real. In what follows, we extend the Bernard-LeClair classification of non-Hermitian matrices [5, 6] to this setting.

If $L$ has a commuting symmetry,

$$LUU^{-1} = L, \qquad (5)$$
then, $L$ and $U$ can be simultaneously diagonalized. Inside each block of $L$, $U$ acts trivially as a multiple of the identity. If no further commuting symmetries exist, we call such a block an irreducible block. The classification of $L$ has to be done for each irreducible block individually. Importantly, the rotation into the diagonal basis of the commuting symmetries $U$ does not in general preserve the Markov properties (2)–(4) in each of the irreducible blocks. For example, most blocks are not represented by a real matrix if $U$ is not itself real with real eigenvalues. Therefore, the classification we develop below is strictly valid only in the cases where no commuting unitary symmetries exist. It is, however, natural to conjecture that the block containing the invariant measure of the Markov process (henceforth referred to as the steady-state block) can still be recast into Markov form. Indeed, consider an initial probability vector fully supported in the steady-state block. Under the Markov dynamics, the initial probability vector evolves into another admissible probability vector. But because all blocks of $L$ are decoupled, no component of the probability vector outside of the steady-state block is generated. We conclude that the steady-state block maps probability vectors into probability vectors and it should be possible to write it in Markov form. In Appendix A we prove this statement under some mild assumptions about the projection into the steady-state block and give an illustrative example of a Markov matrix with SU(2) symmetry. We thus proceed with the classification assuming that either $L$ has no commuting symmetries or, if it has, that we restrict attention to the steady-state block.

We define the (traceless) shifted Markov generator

$$L' = L - \frac{\mathrm{Tr} L}{\mathrm{Tr} 1}, \qquad (6)$$

(it is clear that a commuting symmetry $U$ of $L$ is also one of $L'$ and their irreducible blocks coincide), and look for the possible existence of the following symmetries:

$$SL'S^{-1} = -L', \qquad (7)$$

$$R_+ L'^T R_+^{-1} = L', \qquad (8)$$

$$R_- L'^T R_-^{-1} = -L', \qquad (9)$$

where $S$ and $R_±$ are invertible real matrices only constrained by

$$S^2 = \eta_S, \qquad (10)$$

$$R_+ R_+^{-T} = \eta_+, \qquad (11)$$

$$R_- R_-^{-T} = \eta_-, \qquad (12)$$

$$R_+ R_-^{-T} = \epsilon R_- R_+^{-T}, \qquad (13)$$

where $\eta_S$, $\eta_±$, $\epsilon \in \{\pm 1\}$. Here, $()^T$ denotes transposition, and $R^{-T}$ is a shorthand notation for $R^{-T} := (R^T)^{-1} = (R^{-1})^T$. We do not consider the existence of more than one symmetry of each type, since their composition commutes with $L'$ and we assumed $L'$ to be irreducible. Moreover, the three transformations are not independent, as the composition of any two yields the third. Consequently, we must consider the presence of either only one of them or all three. In the latter case, we take $R_±$ as the independent symmetries and $S = R_+ R_-^{-T}$ as the derived operation.

The constraints (10)–(13) generalize the restrictions imposed on the symmetry operators in quantum mechanics, where $R_±$ and $S$ must be unitary and $S$ always squares to $+1$. On the other hand, we have imposed the reality of $L'$, $R_±$, and $S$, which is absent otherwise. The constraints are derived from the irreducibility of $L'$, by successive application of multiple
<table><caption>**Table 1** Fifteenfold classification of classical Markov matrices</caption><thead><tr><th>Class</th><th>R<sub>+</sub>R<sub>+</sub><sup>-T</sup></th><th>R<sub>-</sub>R<sub>-</sub><sup>-T</sup></th><th>S<sup>2</sup></th><th>Analytical</th><th>Numerical</th><th>Dihedral</th><th>Kramers</th></tr></thead><tbody><tr><td>AI</td><td>0</td><td>0</td><td>0</td><td>Yes</td><td>Yes (f = 0)</td><td>No</td><td>No</td></tr><tr><td>AI<sub>+</sub></td><td>0</td><td>0</td><td>+1</td><td>Yes</td><td>Yes (f = 0)</td><td>Yes</td><td>No</td></tr><tr><td>AI<sub>-</sub></td><td>0</td><td>0</td><td>-1</td><td>No</td><td>Yes (f = 0)</td><td>Yes</td><td>No</td></tr><tr><td>BDI<sup>†</sup></td><td>+1</td><td>0</td><td>0</td><td>Yes</td><td>Yes (f = 0)</td><td>No</td><td>No</td></tr><tr><td>DIII<sup>†</sup></td><td>-1</td><td>0</td><td>0</td><td>No</td><td>Yes (f → 0)</td><td>No</td><td>Yes</td></tr><tr><td>BDI</td><td>0</td><td>+1</td><td>0</td><td>No</td><td>No</td><td>Yes</td><td>No</td></tr><tr><td>CI</td><td>0</td><td>-1</td><td>0</td><td>No</td><td>Yes (f = 0)</td><td>Yes</td><td>No</td></tr><tr><td>BDI<sub>++</sub></td><td>+1</td><td>+1</td><td>+1</td><td>Yes</td><td>No</td><td>No</td><td>No</td></tr><tr><td>CI<sub>+-</sub></td><td>+1</td><td>-1</td><td>+1</td><td>Yes</td><td>Yes (f = 0)</td><td>No</td><td>No</td></tr><tr><td>BDI<sub>-+</sub></td><td>-1</td><td>+1</td><td>+1</td><td>No</td><td>No</td><td>No</td><td>Yes</td></tr><tr><td>CI<sub>++</sub></td><td>-1</td><td>-1</td><td>+1</td><td>No</td><td>Yes (f → 0)</td><td>No</td><td>Yes</td></tr><tr><td>BDI<sub>-+</sub></td><td>+1</td><td>+1</td><td>-1</td><td>No</td><td>No</td><td>No</td><td>No</td></tr><tr><td>CI<sub>--</sub></td><td>+1</td><td>-1</td><td>-1</td><td>No</td><td>Yes (f = 0)</td><td>No</td><td>No</td></tr><tr><td>BDI<sub>-+</sub></td><td>-1</td><td>+1</td><td>-1</td><td>No</td><td>No</td><td>No</td><td>Yes</td></tr><tr><td>CI<sub>-+</sub></td><td>-1</td><td>-1</td><td>-1</td><td>No</td><td>Yes (f → 0)</td><td>No</td><td>Yes</td></tr></tbody></table>

The first column lists the name of the class. (We adapted the nomenclature of Ref. [6] to the present setting; note that while we follow their conventions, our notations differ slightly.) The following three columns give the values of $\eta_+$, $\eta_-$, and $\eta_S$, respectively. (The two indices of the Cartan labels in the class names are, respectively, $\eta_S$ and $\eta_+\eta_-$. ) The column “Analytical” indicates whether we have an analytic parametrization of a family of Markov matrices in the respective class (Sect. 3), while the column “Numerical” indicates whether we have found numerical examples of generators in that class using our Monte Carlo sampling procedure (Sect. 4). For the latter, in parentheses, we indicate whether the cost function of the numerical solution reaches zero in a finite number of steps (“f = 0”) or goes down toward zero with the number of iterations (“f → 0”). The last two columns specify whether the spectra of generators in the respective class have some special properties (Sect. 5), namely, dihedral symmetry (column “Dihedral”) and/or Kramers degeneracy (column “Kramers”)

symmetries, as we now show. To this end, we note that Eqs. (7)–(9) are invariant under $X \mapsto \gamma X$, where $X = S$, $R_\pm$ and $\gamma$ is a real number, and, therefore, we can freely rescale $X$ (we emphasize there is no phase freedom). We start by applying the $S$ symmetry twice: $S^2L'S^2 = L'$. By irreducibility of $L'$, $S^2$ must be a multiple of the identity, $S^2 = \lambda \mathbb{1}$, with $\lambda$ real, and redefining $S \mapsto S/\sqrt{|\lambda|}$, we arrive at Eq.(10). Similarly, acting twice with, say, $R_+$ yields $R_+R_+^T = \lambda \mathbb{1}$ for some real $\lambda$. It follows that $R_+ = \lambda R_+^T = \lambda(\lambda R_+^T)^T = \lambda^2 R_+$ and, hence, $\lambda = \pm 1$, implying Eq. (11). The same applies to $R_-$. Normalizing the “squares” of $R_\pm$ did not require using their scale freedom, which can be used to fix the “commutation relation” of $R_+$ and $R_-$ instead. Indeed, applying them successively to $L'$, we have $R_+(R_-L'T R_-^{-1})^T R_+^{-1} = -L' = R_-(R_+L'T R_+^{-1})^T R_-^{-1}$. The irreducibility of $L'$ then implies that $R_+R_-^T = \lambda R_-R_+^T$ for a real $\lambda$. Redefining $R_+ \mapsto \sqrt{|\lambda|}R_+$, we arrive at Eq. (13). This exhausts the scale freedom of the symmetry operations. Finally, we recall that, in the presence of both $R_\pm$, $S = R_+R_-^T$ is a derived operation. In particular, its square is determined by the “squares” of $R_\pm$ and their “commutation relation”, namely,

$$ \eta_S = \epsilon \eta_+ \eta_-, \qquad (14) $$

and we can thus use them interchangeably.
The symmetry class of (an irreducible block of) $L$ is thus determined by three signed numbers $\pm 1$, the “squares” of the symmetry operations, $S^2$ and $R_{\pm}R_{\pm}^{-T}$. Counting them is now a simple matter (they are summarized in Table 1):

* *No symmetry.* One class, dubbed AI.

* *One symmetry.* We have three types of symmetries, each with two possible values for the “square”, leading to $3 \times 2 = 6$ classes, dubbed $AI_+$, $AI_-$, $BDI^+$, $DIII^+$, $BDI$, and $CI$.

* *Three symmetries.* Each of the three symmetries $R_{\pm}$ and $S$ can “square” to two values, giving $2^3 = 8$ classes, dubbed $BDI_+$, $BDI_-$, $BDI_+$, $BDI_-$, $CI_+$, $CI_-$, $CI_+$, and $CI_-$ (alternatively, as mentioned above, we could specify whether $R_+$ and $R_-$ “commute” or “anticommute” instead of the “square” of $S$).

In total, we thus have $1 + 6 + 8 = 15$ possible symmetry classes in classical stochastic processes. They are in one-to-one correspondence with the 15 non-Hermitian Bernard-LeClair classes with $T_+^2 = -1$ (from which we borrowed their names) [6], but are implemented by nonunitary and nonunitary real matrices that satisfy the generalized constraints (10)–(13).

The allowed set of classes followed simply by restricting our attention to real matrices with real symmetry operations [Eq. (2)] and has not yet taken into account the restrictions of Eqs. (3) and (4). In particular, it is not guaranteed that all 15 classes can be realized with those restrictions. In the next two sections, we look for explicit examples of Markov matrices in each of these classes. First, in Sect. 3 we will follow the usual procedure [9, 14, 15] and look for analytic parameterizations of families of matrices of arbitrary dimension. However, because of the Eqs. (3) and (4), we can only do so in five of the fifteen classes. Therefore, we then turn to obtain further examples by Monte Carlo sampling in Sect. 4, which yields examples in ten classes. Remarkably, there is one class with an analytic expression that is *not* found by our numerical sampling algorithm. Therefore, we establish that the symmetry classification of Markov generators subjected to Eqs. (3) and (4) is *at least* elevenfold but because of this limitation of the Monte Carlo algorithm, we cannot definitively rule out the existence of the remaining four.

## 3 Analytical Families of Markov Matrices

We first give explicit realizations of Markov matrices for five classes, namely, $AI$, $AI_+$, $BDI^+$, $BDI_+$, and $CI_+$. These simple descriptions are possible because in those cases there is no mismatch between the preferred computational basis, in which the positivity constraints are simple, and the basis in which the operators implementing the discrete symmetries are simple. Borrowing the notation to be introduced in Sect. 4, these cases correspond to those for which $W = 1$.

To this end, we parametrize the classical Markov generator $L$ in terms of a real matrix $M$ with arbitrary diagonal entries and nonnegative off-diagonal entries:

$$L_{ij} = M_{ij} - \delta_{ij} \sum_{k=1}^{N} M_{kj}, \quad M_{ij} \ge 0, \quad i \neq j. \qquad (15)$$
Without loss of generality, we take M to be traceless. To ensure that L has the appropriate symmetry, we demand that M satisfies the S or $R_±$ symmetries, accordingly:

$$SMS^{-1} = -M, \qquad (16)$$

$$R_{-}M^{T}R_{-}^{-1} = -M, \qquad (17)$$

$$R_{+}M^{T}R_{+}^{-1} = +M. \qquad (18)$$

However, this is not sufficient to guarantee a symmetry of L. We further impose the following necessary condition:

$$\sum_{n} M_{nk} = -\frac{\alpha}{2}, \quad \forall k, \qquad (19)$$

with $\alpha = \frac{2}{N} \text{Tr} L$. Equation (19) is a very restrictive condition as it demands that summing matrix elements of M on each column yields the same output. Notwithstanding, it amounts only to a gauge choice ensuring that a symmetry of M is passed to L and does not restrict the physical processes. To see this, notice that we can shift M by any diagonal matrix, i.e., $M_{ij} \rightarrow \tilde{M}_{ij} = M_{ij} - \delta_{ij}x_j$ gives the same matrix L. The gauge degrees of freedom $x_j$ are exactly enough to always impose that both $\text{Tr} M = 0$ and $\sum_k M_{kj} = -\alpha/2$. Finally, note that, assuming Eq. (19), the constant $\alpha$ is determined by summing over $k$: since $\text{Tr} M = 0$ by gauge choice, Eq. (15) implies $\text{Tr} L = -\sum_{kj} M_{kj}$.

Let us first consider the case of an S symmetry. We have:

$$\begin{aligned} \sum_{kl} S_{ik} L_{kl} (S^{-1})_{lj} &= \sum_{kl} S_{ik} M_{kl} (S^{-1})_{lj} - \sum_{kl} S_{ik} \delta_{kl} \sum_{n} M_{nl} (S^{-1})_{lj} \\ &= -M_{ij} - \sum_{k} S_{ik} (S^{-1})_{kj} \sum_{n} M_{nk}. \end{aligned} \qquad (20)$$

To have an S symmetry, we need, in addition to Eq. (16), that

$$\sum_{k} S_{ik} (S^{-1})_{kj} \sum_{n} M_{nk} = -\delta_{ij} \sum_{n} M_{nj}, \qquad (21)$$

which is achieved if Eq. (19) holds. Then, Eq. (20) reads as:

$$\sum_{kl} S_{ik} L_{kl} (S^{-1})_{lj} = -M_{ij} + \frac{\alpha}{2} \delta_{ij} = -M_{ij} + \delta_{ij} \sum_{n} M_{nj} + \alpha \delta_{ij} = -L_{ij} + \alpha \delta_{ij}, \qquad (22)$$

which is Eq. (7) for the shifted Markov matrix.

Next, we consider the case of a $R_±$ symmetry. As before, we have:

$$\sum_{kl} (R_{\pm})_{ik} L_{kl} (R_{\pm}^{-1})_{lj} = \pm M_{ij} - \sum_{k} (R_{\pm})_{ik} (R_{\pm}^{-1})_{kj} \sum_{n} M_{nk}. \qquad (23)$$

We notice that the same condition (19) is sufficient to guarantee the $R_±$ symmetry: for $R_-$, the reasoning is exactly the same as above; for $R_+$, we note that

$$\sum_{kl} (R_{+})_{ik} L_{kl} (R_{-}^{-1})_{lj} = M_{ij} + \frac{\alpha}{2} \delta_{ij} = M_{ij} - \delta_{ij} \sum_{n} M_{nj} = L_{ij}. \qquad (24)$$

We can now follow the procedure of Ref. [15] and seek the simplest irreducible parametrization of the elements of each symmetry class that is consistent with the symmetries of that class, Eqs. (16)–(18), and the positivity of the entries. By irreducible we again mean that no additional unitary symmetries exist. Equation (19) has to be imposed on top
of the parametrizations we derive. Throughout, X, Y, Z denote the three real traceless 2 × 2 matrices

$$X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad Y = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}, \quad Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \qquad (25)$$

1 is the N-dimensional identity and A, B, C, ... are N × N real matrices (arbitrary, unless specified otherwise).

### 3.1 Class AI, $S^2 = R_+R_+^{-T} = R_-R_-^{-T} = 0$

There are no symmetries to implement. M is an arbitrary real matrix with nonnegative off-diagonal entries.

### 3.2 Class AI+, $S^2 = +1, R_+R_+^{-T} = R_-R_-^{-T} = 0$

Class AI+ is labeled by the symmetry $S^2 = +1$. The simplest nontrivial choice of S satisfying $S^2 = +1$ and compatible with the reality of the Markov matrix is $S = Z \otimes 1$. In this case, M is an $2N \times 2N$ matrix with block structure

$$M = \begin{pmatrix} A & B \\ C & D \end{pmatrix}. \qquad (26)$$

Equation (16) reads as

$$SMS^{-1} = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} = \begin{pmatrix} A & -B \\ -C & D \end{pmatrix} = \begin{pmatrix} -A & -B \\ -C & -D \end{pmatrix} = -M, \qquad (27)$$

from which it follows that the diagonal blocks vanish ($A = D = 0$) because of the positivity of the entries, while there are no further conditions on the off-diagonal blocks (and we can rename B to A and C to B).

To apply Eq. (19), we define

$$\tilde{J}_{ij} = \frac{J_{ij}}{\sum_k J_{kj}} \qquad (28)$$

for an arbitrary matrix J. The parametrization of a Markov generator in class AI+ is finally found to be:

$$M = \begin{pmatrix} 0 & \tilde{A} \\ \tilde{B} & 0 \end{pmatrix}. \qquad (29)$$

### 3.3 Class BDI$^†$, $R_+R_+^{-T} = +1, R_-R_-^{-T} = S^2 = 0$

Class BDI$^†$ is labeled by the symmetry $R_+R_+^{-T} = +1$. Taking $R_+ = X \otimes 1$ and M of the form (26), Eq. (18) reads as

$$R_+ M^T R_+^{-1} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} A^T & C^T \\ B^T & D^T \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} D^T & B^T \\ C^T & A^T \end{pmatrix} = \begin{pmatrix} A & B \\ C & D \end{pmatrix} = M, \qquad (30)$$

whence it follows that

$$M = \begin{pmatrix} A & B \\ C & A^T \end{pmatrix}, \quad \left\{ \begin{array}{l} B = B^T \\ C = C^T \end{array} \right. . \qquad (31)$$
Now, to apply Eq. (19), we cannot simply renormalize the matrix elements by replacing A, B, and C by $\tilde{A}$, $\tilde{B}$, and $\tilde{C}$ because the transformation (28) spoils the transposition symmetry ($\tilde{B} \neq \tilde{B}^T$, etc.). Instead, we can proceed as follows. We can choose the entries of the upper-left block A freely. This also fixes the lower-right block $A^T$. We can also choose arbitrary upper-triangular entries for B and C; $B = B^T$ and $C = C^T$ then fixes the lower-triangular entries. Finally, the diagonal entries of B and C are chosen such that all columns sum to the same value (this fixes all but one diagonal element, which can be chosen freely).

### 3.4 Class $CI_{+-}$, $R_+R_+^T = S^2 = +1$, $R_-R_-^T = -1$

Class $CI_{+-}$ is labeled by the symmetries $R_+R_+^T = S^2 = +1$ and $R_-R_-^T = -1$. We can take the parametrization of $AI_+$ matrix, which has the $S^2 = +1$ symmetry and, on top, require the $R_-R_-^T = -1$ symmetry. We implement the $R_-$ symmetry as $R_- = Y \otimes 1$. Equation (17) reads as

$$R_- M^T R_-^{-1} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} 0 & B^T \\ A^T & 0 \end{pmatrix} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} 0 & -A^T \\ -B^T & 0 \end{pmatrix} = \begin{pmatrix} 0 & -A \\ -B & 0 \end{pmatrix} = -M, \quad (32)$$

which implies that $A = A^T$ and $B = B^T$. Finally, the existence of the $S$ and $R_-$ symmetries automatically implies the existence of the third one, $R_+$, which is given by $R_+ = X \otimes 1$. Indeed, the matrix $M$ we found above satisfies Eq. (18):

$$R_+ M^T R_+^{-1} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 0 & B^T \\ A^T & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} 0 & A^T \\ B^T & 0 \end{pmatrix} = \begin{pmatrix} 0 & A \\ B & 0 \end{pmatrix} = M. \quad (33)$$

Note that this corresponds to Eq. (31) with $A = 0$ (because we need vanishing diagonal blocks to implement a $S^2 = 1$ symmetry). Hence, to apply Eq. (19), we can follow the same procedure as described there. The parameterization of a Markov generator in class $CI_{+-}$ is thus:

$$M = \begin{pmatrix} 0 & A \\ B & 0 \end{pmatrix}, \quad \begin{cases} A = A^T \\ B = B^T \end{cases}, \quad (34)$$

where the diagonal elements of A and B are fixed as described after Eq. (31).

### 3.5 Class $BDI_{++}$, $R_+R_-^T = R_-R_-^T = S^2 = +1$

The final example is that of class $BDI_{++}$, labeled by the symmetries $R_+R_+^T = R_-R_-^T = S^2 = +1$, which requires a 4 × 4 block structure. We take $S = Z \otimes Z \otimes 1$, $R_+ = X \otimes X \otimes 1$, and $R_- = Y \otimes Y \otimes 1$. Following the same procedure as before, we arrive at

$$M = \begin{pmatrix} 0 & A & C & 0 \\ B & 0 & 0 & C^T \\ D & 0 & 0 & A^T \\ 0 & D^T & B^T & 0 \end{pmatrix}. \quad (35)$$

Finally, to apply Eq. (19), we note that it can be written as a linear system of $4N$ equations (of which only $4N-1$ are independent) with $4N^2$ entries of A, B, C, and D as variables. We take the dependent variables to be $B_{NN}$, $C_{iN}$, $C_{Ni}$, $D_{iN}$, and $D_{Ni}$, $i = 1, \dots, N$. Choosing a large enough positive value for the sum of the columns, we can fix all other entries freely, and then solve for the dependent variables.
## 3.6 Absence of Analytic Examples for the Other Classes

For the remaining classes, the positivity constraint sets too many entries of $M$ to zero, which induces additional symmetries: either (i) originally-absent $S$ and $R_±$ symmetries, or (ii) commuting symmetries that break the desired $S$ or $R_±$ symmetry; in either case, the symmetry class is changed.

As an illustrative example of case (i), we consider class CI, characterized by $R_-R_-^{-T} = -1$, $R_+R_+^{-T} = S^2 = 0$. If we implement $R_- = Y \otimes 1$, then, proceeding as before, we find the matrix representation

$$ M = \begin{pmatrix} A & B \\ C & -A^T \end{pmatrix}, \quad \begin{cases} B = B^T \\ C = C^T \end{cases}. \qquad (36) $$

Because of the positivity of the entries, we must set $A = 0$. However, this reduces Eqs. (36) to (34), i.e., induces additional symmetries $R_+ = X \otimes 1$ and $S = Z \otimes 1$, and gives, instead a parametrization of class $AI_{+-}$.

As an example of case (ii), we turn to class $AI_-$, characterized by $S^2 = -1$ and $R_+R_+^{-T} = R_-R_-^{-T} = 0$. Choosing $S = Y \otimes 1$, we obtain

$$ M = \begin{pmatrix} 0 & A \\ A & 0 \end{pmatrix}. \qquad (37) $$

While this matrix has no additional $R_±$ symmetries (which would change its class), it admits a second $S$ symmetry, given by $S' = Z \otimes 1$ (as is evident from the chiral structure of $M$). Now, as discussed in Sect. 2, two symmetries of the same type lead to a matrix that commutes with $M$, this case $U = X \otimes 1$ (which, up to a phase, gives the product $SS'$, and, consequently, $M$ block diagonalizes in the basis where $U$ is diagonal). Crucially, $U$ *anticomutes* with both $S$ and $S'$ and, hence, $S$ and $S'$ connect states with opposite eigenvalue of $U$ and do not act inside a single block of $M$. As a consequence, the irreducible blocks of $L$ do not have any symmetry and belong to class $AI$ instead of $AI_-$.

The discussion of the preceding two examples can be systematized. We have checked that for operators $R_±$ and $S$ implemented by any string of up to four Pauli matrices (corresponding to matrices with up to a 16 × 16 block structure), no classes beyond the five discussed in Sects. 3.1–3.5 can be built in the computational basis.

# 4 Systematic Numerical Solutions of the Symmetry Classes

We have seen that, working in the computational basis, only five of the fifteen classes of Markov matrices can be realized. In the following, we look for the existence of a different basis where $L$ has real and positive entries. To this end, we develop a simple but systematic stochastic optimization algorithm.

## 4.1 Parametrization of the Symmetry Operators

Let us start by choosing a convenient general parametrization of the symmetry operators $S$ and $R_±$. We first consider the classes with $S$ symmetry only. Without loss of generality, because $S^2 = \pm 1$, we can write $S = W\Sigma_S W^{-1}$, where $W$ is an arbitrary invertible matrix
<table><caption>**Table 2** Parametrization of the symmetry operators for the eight classes with all three symmetries</caption><thead><tr><th>Class</th><th>η<sub>+</sub></th><th>η<sub>-</sub></th><th>ε</th><th>Σ<sub>+</sub></th><th>Σ̃<sub>-</sub></th></tr></thead><tbody><tr><td>BDI<sub>++</sub></td><td>+1</td><td>+1</td><td>+1</td><td>diag(s<sub>1</sub>,...,s<sub>N</sub>)</td><td>diag(s<sub>1</sub>,...,s<sub>N</sub>)</td></tr><tr><td>CI<sub>-</sub></td><td>+1</td><td>-1</td><td>+1</td><td>s<sub>1</sub>I ⊕ ... ⊕ s<sub>N/2</sub>I</td><td>s<sub>1</sub>Y ⊕ ... ⊕ s<sub>N/2</sub>Y</td></tr><tr><td>BDI<sub>-</sub></td><td>-1</td><td>+1</td><td>+1</td><td>s<sub>1</sub>Y ⊕ ... ⊕ s<sub>N/2</sub>Y</td><td>s<sub>1</sub>I ⊕ ... ⊕ s<sub>N/2</sub>I</td></tr><tr><td>CI<sub>++</sub></td><td>-1</td><td>-1</td><td>+1</td><td>s<sub>1</sub>Y ⊕ ... ⊕ s<sub>N/2</sub>Y</td><td>s<sub>1</sub>Y ⊕ ... ⊕ s<sub>N/2</sub>Y</td></tr><tr><td>BDI<sub>-</sub></td><td>+1</td><td>+1</td><td>-1</td><td>s<sub>1</sub>Z ⊕ ... ⊕ s<sub>N/2</sub>Z</td><td>s<sub>1</sub>X ⊕ ... ⊕ s<sub>N/2</sub>X</td></tr><tr><td>CI<sub>+-</sub></td><td>+1</td><td>-1</td><td>-1</td><td>s<sub>1</sub>Z ⊕ ... ⊕ s<sub>N/2</sub>Z</td><td>s<sub>1</sub>Y ⊕ ... ⊕ s<sub>N/2</sub>Y</td></tr><tr><td>BDI<sub>+-</sub></td><td>-1</td><td>+1</td><td>-1</td><td>s<sub>1</sub>Y ⊕ ... ⊕ s<sub>N/2</sub>Y</td><td>s<sub>1</sub>Z ⊕ ... ⊕ s<sub>N/2</sub>Z</td></tr><tr><td>CI<sub>-</sub></td><td>-1</td><td>-1</td><td>-1</td><td>s<sub>1</sub>(Y ⊗ I) ⊕ ... ⊕ s<sub>N/4</sub>(Y ⊗ I)</td><td>s<sub>1</sub>(Z ⊗ Y) ⊕ ... ⊕ s<sub>N/4</sub>(Z ⊗ Y)</td></tr></tbody></table>

with fixed norm, Tr $W^T W = 1$, and

$$ \left\{ \begin{array}{ll} \Sigma_S = \operatorname{diag}(+1, \dots, +1, \underbrace{-1, \dots, -1}_{n}, \underbrace{-1, \dots, -1}_{N-n}), & \text{if } \eta_S = +1, \\ \Sigma_S = \underbrace{Y \oplus \dots \oplus Y}_{N/2}, & \text{if } \eta_S = -1, \end{array} \right. \qquad (38) $$

with diag denoting a diagonal matrix.

If, instead, there is only a $R_\pm$ symmetry, implemented by either a real symmetric ($\eta_\pm = +1$) or a real antisymmetric ($\eta_\pm = -1$) matrix, without loss of generality, we can write it as $R_\pm = W\Sigma_\pm W^T$, with W now orthogonal and

$$ \left\{ \begin{array}{ll} \Sigma_\pm = \operatorname{diag}(s_1, \dots, s_N) & \text{if } \eta_\pm = +1, \\ \Sigma_\pm = s_1 Y \oplus \dots \oplus s_{N/2} Y & \text{if } \eta_\pm = -1, \end{array} \right. \qquad (39) $$

where $s_j$ are real numbers. Because of the homogeneity of the symmetry conditions, we can fix Tr $R_\pm R_\pm^T = \sum_j s_j^2 = 1$.

Finally, we turn to the case with all three symmetries. We have seen that, in that case, we do not need to consider $S$ explicitly, as its square is contained in the compatibility condition (13). In this case, we parametrize $\Sigma_+$ as in Eq. (39) and define $\Sigma_-$ through $R_- = W\Sigma_- W^T$. Note that, because we are using the same $W$, $\Sigma_-$ does not need to satisfy Eq. (39), but it still satisfies $\Sigma_-^T = \eta_-\Sigma_-$. Plugging the definitions of $\Sigma_\pm$ into Eq. (13) yields

$$ \Sigma_+ \Sigma_-^{-T} = \epsilon \Sigma_- \Sigma_+^{-T}. \qquad (40) $$

For example, for class BDI<sub>+</sub> ($R_+R_+^T = R_-R_-^T = +1$, $S^2 = -1$), Eq. (40) requires that only $N/2$ $s_j$ are independent, and that $\Sigma_\pm$ take the form $\Sigma_+ = s_1Z \oplus \dots \oplus s_{N/2}Z$ and $\Sigma_- = g\tilde{\Sigma}_-$ with $\tilde{\Sigma}_- = s_1X \oplus \dots \oplus s_{N/2}X$. Here, $g = (\pm I) \oplus \dots \oplus (\pm I)$, where $I$ the two-dimensional identity and all signs can be independently chosen. By reorganizing basis elements, $g$ can be brought to the form $g = \operatorname{diag}(+1, \dots, +1, \underbrace{-1, \dots, -1}_{n}, \underbrace{-1, \dots, -1}_{N/2-n}) \otimes I$, so that we

only need to specify the number $n$ of $+1$ entries. We have also assumed that all independent $s_j$ are nondegenerate. We can proceed similarly to find the general parametrization of the symmetry operators in the remaining classes. The result for the eight classes with all three symmetries is listed in Table 2.
## 4.2 Solution of the Symmetry Constraints

For fixed symmetry operators $S$ and $R_±$, the symmetry transformations (7)–(9) together with the Markov condition (4) form a linear system for the entries of $L$ that we can solve exactly. To this end, we (row-) vectorize

$$L \mapsto |L\rangle = (L_{11}, L_{12}, \dots, L_{1N}, L_{21}, \dots, L_{NN})^T, \quad (41)$$

in which case, Eqs. (7)–(9) read as, respectively,

$$(S \otimes S^{-T} + \mathbb{1}_{N^2} - \frac{2}{N} |\mathbb{1}_N\rangle\langle\mathbb{1}_N|) |L\rangle = 0, \quad (42)$$

$$((R_+ \otimes R_+^{-T})S - \mathbb{1}_{N^2}) |L\rangle = 0, \quad (43)$$

$$((R_- \otimes R_-^{-T})S + \mathbb{1}_{N^2} - \frac{2}{N} |\mathbb{1}_N\rangle\langle\mathbb{1}_N|) |L\rangle = 0, \quad (44)$$

where $S$ is the SWAP operator, $S(A \otimes B)S = B \otimes A$ for arbitrary operators $A$ and $B$, which implements transposition of $L$. On the other hand, Eq. (4) becomes

$$\left( \mathbf{e}_N^T \otimes \mathbb{1}_N \right) |L\rangle = 0, \quad (45)$$

where $\mathbf{e}_N$ is the $N$-dimensional (column) vector with all entries equal to one. We can write the combined system of equations as $\mathcal{F}|L\rangle = 0$. This is an overdetermined homogeneous system of equations with $N^2$ variables and $N^2 + N$ (for classes with one symmetry) or $N^2 + 2N$ (for classes with three symmetries) equations (recall that, in the case of three symmetries, Eq. (42) is already implied by Eqs. (43) and (44) and can be omitted). Nevertheless, for large enough $N$ ($N \ge 2, 4$, or $6$ depending on the class), there is always a nonzero family of solutions. However, a generic member of this family does not satisfy the positivity condition (3). It is not even guaranteed that, for a given $\mathcal{F}$ (i.e., given $S$ and $R_±$), there exists any member satisfying it.

To proceed, we construct a basis $\{|b_j\rangle\}_{j=1}^K$ for the $K$-dimensional kernel of $\mathcal{F}$. We write a solution of the homogeneous system as $|L\rangle = \sum_{j=1}^K \ell_j |b_j\rangle$, where, without loss of generality, we can fix $\sum_j \ell_j^2 = 1$. We introduce a cost function

$$f := \frac{1}{N} \sum_{j \neq k} \frac{|L_{jk}| - L_{jk}}{2}, \quad (46)$$

which satisfies $f = 0$ for a Markov matrix and $f > 0$ for a non-Markov matrix, with a larger value signaling a larger departure from Markovianity, and look for the $\ell$ that minimizes it:

$$\ell_* = \underset{\ell \in S^K}{\operatorname{argmin}} f(\ell). \quad (47)$$

## 4.3 Monte Carlo Sampling

We now use a stochastic optimization algorithm to minimize the cost function $f$. For a fixed class, we start from a random initial set $\{W^{(0)}, s^{(0)}\}$ and a choice of $g$. We use a (zero-temperature) Metropolis algorithm to decrease $f$, where at each step we perturb $\{W, s\}$ with
strength $\delta$. More precisely, we perturb

$$s^{(i)} \rightarrow s^{(i+1)} = \frac{s^{(i)} + \delta v}{|s^{(i)} + \delta v|}, \qquad (48)$$

where $v$ is a vector of Gaussian random variables with zero mean and unit variance. $W$ is perturbed differently depending on whether it is orthogonal (classes with an $R_\pm$ symmetry) or not (classes with only $S$ symmetry). In the former, we use

$$W^{(i)} \rightarrow W^{(i+1)} = e^{\delta A} W^{(i)}, \qquad (49)$$

where $A$ is a Gaussian real antisymmetric random matrix, whereas in the latter, we set

$$W^{(i)} \rightarrow W^{(i+1)} = \frac{W^{(i)} + \delta V}{||W^{(i)} + \delta V||}, \qquad (50)$$

where $V$ is a real Gaussian random matrix and $|| \cdot ||$ is the Hilbert-Schmidt norm. For each $\{W^{(i)}, s^{(i)}\}$, we compute the cost function $f$ and, if it is lower than that of the previous iteration, we accept the new configuration. After $m_\delta$ iterations without an acceptance, we decrease $\delta \rightarrow \delta/2$.

## 4.4 Results

Fig. 1 depicts the results of the stochastic optimization algorithm for $N = 8$, $n = 2$ for all classes except AI (which has no symmetry constraints and therefore no algorithm to run). We started the algorithm with $\delta = 1$, ran it for $2 \times 10^4$ steps (or until $f < 10^{-12}$ was reached), and chose $m_\delta = 500$. In nine of the classes (plus AI), there are solutions that satisfy the positivity condition. Of these, for six—AI+, AI-, BDI†, CI, CI+, and CI——$f = 0$ is reached exactly after a finite (and small) number of iterations of the algorithm. For the other three— DIII†, CI++, and CI+—$f$ decreases arbitrarily with the number of iterations but we have not reached $f < 10^{-12}$ for any of our runs. Nevertheless, the results for these three classes are strikingly different from the five classes—BDI, BDI++, BDI+, BDI+, and BDI—— where $f$ plateaus at a constant of order 1 and does not decrease with the number of iterations. We interpret this as a solution that is reached asymptotically for the former three classes, while no solution is found for the latter five. Remarkably, for one of the classes where we could not find a numerical solution (namely, BDI++), there is an analytical solution given in Sect. 3.5. This shows that not finding a numerical solution is a limitation of the algorithm and does not strictly rule out the existence of the four remaining classes. We checked extensively that different choices of sampling of $s$ and $W$, values of $n$, acceptance criteria (including finite-temperature algorithms with different annealing schedules and varying initial values of $\delta$ and $m_\delta$), cost functions, and allowing for complex entries lead to qualitatively similar results; in particular, the classes for which solutions exist or not and whether $f = 0$ is reached asymptotically or in a finite number of steps did not depend on any of these choices.

We conclude with some empirical observations extracted from the results in Fig. 1 and Table 1: no solution was found in any class with $R_-R_-^T = +1$; in all classes with $R_-R_-^T \neq +1$ and $R_+R_+^T = -1$, the solution was reached asymptotically; in all other classes (i.e., those with $R_-R_-^T \neq +1$ and $R_+R_+^T \neq -1$), a solution was found exactly after a finite number of iterations. While we do not have a theoretical understanding of these features of the algorithm, they could help in devising new ones that could, potentially, show the existence of examples in the four missing classes.
<figure><img src="image_2.png" /><figcaption>Fig. 1 Results of the stochastic optimization algorithm for all fourteen classes with symmetry constraints, for $N = 8$ and $n = 2$. Each panel shows the value of the cost function $f$ as a function of the number of acceptances starting from 10 different initial conditions {$W^{(i)}, s^{(i)}$}, with $\delta = 1, m_{\delta} = 500$, and up to $2 \times 10^4$ iterations</figcaption></figure>

# 5 Examples and Physical Consequences

## 5.1 Example: Bipartite Graph

We now discuss simple physical examples of the families built in Sect. 3. The simplest example of a Markov generator in a nontrivial class is that of an exclusion process on a bipartite graph. Let us consider a graph whose $2N$ nodes are either black ($\bullet$) or white ($\circ$) and that there is no hopping between two nodes of the same type. Let us also denote the $j$th node of type $\bullet$ ($\circ$) as $\bullet_j$ ($\circ_j$). In the absence of any further constraints, $S = Z \otimes 1$ plays the role of sublattice symmetry and the Markov matrix of the stochastic process is given by

$$M = \begin{pmatrix} 0 & A \\ B & 0 \end{pmatrix}, \qquad (51)$$
where A and B are $N \times N$ real matrices containing the transition rates $\bullet \to \circ$ and $\circ \to \bullet$, respectively (the transitions $\bullet \to \circ$ and $\circ \to \bullet$ are not necessarily symmetric). The sum of the $k$th column of $M$ gives the rate at which particles hop into the $k$th node of the graph, coming from any other node. Imposing Eq. (19) thus amounts to fixing this total incoming rate into a node to be the same for all nodes. The generator $L$ then has $S^2 = +1$ symmetry. In the absence of any further restrictions, $L$ belongs to class $AI_+$.

If, we further impose that the hopping rate $\bullet_j \to \circ_k$ is the same as $\bullet_k \to \circ_j$, and similarly, the rates $\circ_j \to \bullet_k$ and $\circ_k \to \bullet_j$ coincide, then $A = A^T$ and $B = B^T$. (Note that we are not assuming any symmetry between the hoppings $\bullet_j \to \circ_k$ and $\circ_j \to \bullet_k$, i.e., $A \neq B^T$.) Thus, in addition to the $S$ symmetry, we also have $R_+ = X \otimes 1$, satisfying $R_+ R_+^T = +1$. The simultaneous existence of $S$ and $R_+$ symmetries automatically implies the existence of a $R_-$ symmetry, implemented here by $R_- = Y \otimes 1$, $R_- R_-^T = -1$. Keeping the total incoming rate into any node still fixed, we have a Markov generator of class $CI_{++}$.

If we allow for hopping between nodes of the same type, but which are restricted by the special property that the rate $\bullet_j \to \circ_k$ is equal to the (reverse) rate $\circ_k \to \circ_j$, then we can write the Markov matrix $M$ as

$$ M = \begin{pmatrix} C & A \\ B & C^T \end{pmatrix}, \quad \left\{ \begin{array}{l} A = A^T \\ B = B^T \end{array} \right. , \qquad (52) $$

where $C$ contains the $\bullet \to \circ$ and $\circ \to \bullet$ rates. The diagonal blocks break the $S$ and $R_-$ symmetries, but preserve the $R_+$ symmetry, provided that the total incoming rate into every node is still the same. In that case, the Markov generator $L$ belongs to class $BDI^\dagger$.

Finally, we would need to consider a quadripartite graph to get the other class with an explicit parametrization, $BDI_{++}$.

## 5.2 Dihedral Symmetry of the Spectrum

Having presented some examples of stochastic processes in different classes in the previous section, we now turn to the spectral consequences of different symmetries. Because $L$ is always real, its eigenvalues are real or come in complex conjugate pairs. Likewise, the presence of an $S$ or $R_\pm$ symmetry also pairs different eigenvalues, as we now show.

Let $\lambda_j, j = 1, \dots, N$, be the eigenvalues of $L'$ and $\phi_j$ and $\tilde{\phi}_j$ the corresponding right (column) and left (row) eigenvectors, i.e.,

$$ L'\phi_j = \lambda_j\phi_j, \qquad (53) $$

$$ \tilde{\phi}_j L' = \lambda_j \tilde{\phi}_j \qquad (54) $$

Since $L'$ is real, all eigenvalues are either real or come in pairs $(\lambda_j, \lambda_j^*)$. Likewise, multiplying Eq. (53) by $S$ from the left and using Eq. (7), we have

$$ L'S\phi_j = -\lambda_j S\phi_j, \qquad (55) $$

that is, $S$ reflects the spectrum across the origin since $S\phi_j$ is an eigenvector of $L'$ with the eigenvalue $-\lambda_j$. Consequently, all eigenvalues are zero or come in pairs $(\lambda_j, -\lambda_j)$. Combined with the reflection of the spectrum across the real axis due to the reality of $L'$, it implies dihedral symmetry of the spectrum of $L'$, analogous to that observed in the quantum case [9, 16, 17].
<figure><img src="image_3.png" /><figcaption>Fig. 2 Spectrum of L in the complex plane for a random Markov matrix from class $AI_+$ (left), sampled according to Eq. (29) and which displays dihedral symmetry, and class $BDI^\dagger$ (right), sampled from Eq. (31) and which does not show dihedral symmetry. In both cases, we set $N = 100$ and the entries of the blocks A, B, C are the absolute values of normal random variables with zero mean and variance $2/N$.</figcaption></figure>

On the other hand, multiplying Eq. (53) by $R_\pm$ from the left and using Eqs. (8) or (9), we have

$$ (R_\pm \phi_j)^T L' = \pm \lambda_j (R_\pm \phi_j)^T, \quad (56) $$

i.e., $(R_\pm \phi_j)^T$ is a left eigenvector corresponding to the eigenvalue $\pm \lambda_j$. For $R_-$, eigenvalues thus always come in pairs $(\lambda_j, -\lambda_j)$, and as in the case of an S symmetry, the spectrum displays dihedral symmetry. For $R_+$ no dihedral symmetry is present, as it relates eigenvectors with the same eigenvalues.

To illustrate these properties, we show in Fig. 2 the spectra of a representative member of two classes with an analytical parametrization (see Sect. 3), namely, one with dihedral symmetry (class $AI_+$) and one without (class $BDI^\dagger$).

## 5.3 Kramers Degeneracy

Combining the results of the previous section about the reality of $L'$ and its transformation under $R_+$, it follows that for a given right eigenvector $\phi_j$, $(R_+ \phi_j)^T$ is a left eigenvector corresponding to the same eigenvalue $\lambda_j$. But $\tilde{\phi}_j$ is also a left eigenvector associated with the eigenvalue $\lambda_j$. Now, by biorthogonality, we have $\tilde{\phi}_j \phi_j = 1$. On the other hand,

$$ (R_+ \phi_j)^T \phi_j = (R_+ \phi_j)^T (R_+^{-1} R_+ \phi_j) = (R_+^{-T} R_+ \phi_j)^T (R_+ \phi_j) = \eta_+ \phi_j^T (R_+ \phi_j), \quad (57) $$

whence

$$ (R_+ \phi_j)^T \phi_j = 0, \quad \text{if } \eta_+ = -1. \quad (58) $$

Therefore, $\tilde{\phi}_j$ and $(R_+ \phi_j)^T$ are linearly independent and every eigenvalue $\lambda_j$ is doubly degenerate, a result known as Kramers degeneracy in quantum theory and which we have shown to extend to the classical stochastic setting even when the antiunitarity conditions are relaxed.

None of the examples for which we have analytical parametrization have $R_+ R_+^T = -1$. Nevertheless, the numerical examples written out explicitly in Appendices B.2, B.4, and B.6 show double degeneracy of all of their eigenvalues, as expected.
## 5.4 Detailed Balance and $R_+$ Symmetry

A Markov matrix $L$ satisfies detailed balance if there exist positive numbers $\pi_j$ such that

$$L_{ij} \pi_j = L_{ji} \pi_i. \qquad (59)$$

Since this equation is linear in $\pi$, we can assume without loss of generality that $\sum_j \pi_j = 1$. If Eq. (59) holds, $\pi$ forms the stationary measure of the process defined by $L$ (because $\partial_t \pi_i = \sum_j L_{ij} \pi_j = \pi_i \sum_j L_{ji} = 0$). Clearly, a Markov matrix satisfying the detailed balance of Eq. (59) possesses an $R_+$ symmetry as it is enough to choose

$$R_+ := \operatorname{diag}(\pi_1, \pi_2, \dots, \pi_N) \qquad (60)$$

to ensure Eq. (8), i.e. $R_+ L^T = L R_+$ (recall that for $R_+$ the shift in $L$ is irrelevant so that Eq. (8) can be written in terms of $L$ instead of $L'$).

More generally, we can define a vector $p_* := R_+ e_N$ (where, as before, $e_N$ is the $N$-dimensional vector with all entries equal to 1), with component $(p_*)_i = \sum_j (R_+)_{ij}$. By construction, $p_*$ is an invariant state of $L$ (because $\partial_t p_* = LR_+e_N = R_+L^T e_N = 0$). Assuming that $p_*$ is the unique invariant state of $L$ (which should be the case for a generic Markov process), then, by the Perron-Frobenius theorem, all its components are positive, and $p_*$ is an invariant measure. Thus, an $R_+$ symmetry can be thought of as nondiagonal detailed balance.

Finally, we note that detailed balance, Eq. (59) can be seen as a special case of self-duality of the Markov process described by $L$ [18–20]. More generally, duality of Markov processes can be used, e.g., to write $m$-point correlation functions in $N$-particle systems in terms of correlation functions in $m$-particle systems [18, 19].

## 5.5 Time Reversal of Stochastic Trajectories

Consider a Markov chain with states $i, j, \dots$ and the stochastic matrix $P = e^{\delta t L}$, for some time interval $\delta t$, whose entries are the transition probabilities $P_{ij} = \operatorname{Prob}(i \to j)$. The probability to follow a given stochastic trajectory $\omega = (i_0, i_1, \dots, i_T)$ or its reverse $\omega \leftarrow = (i_T, i_{T-1}, \dots, i_0)$ are given by, respectively,

$$\mathbb{P}(\omega) = P_{i_0 i_1} P_{i_1 i_2} \cdots P_{i_{T-1} i_T} \quad \text{and} \quad \mathbb{P}(\omega \leftarrow) = P_{i_1 i_0} P_{i_2 i_1} \cdots P_{i_T i_{T-1}}. \qquad (61)$$

In the presence of detailed balance, Eq. (59), which is inherited by $P$, we can immediately relate the trajectory $\omega$ to its reverse through:

$$\mathbb{P}(\omega \leftarrow) = \frac{\pi_{i_T}}{\pi_{i_0}} \mathbb{P}(\omega). \qquad (62)$$

This remark can simply be extended to the cases in which $R_+$ is an involutive permutation, up to a redefinition of the time-reversed trajectory.

In the preceding discussion, we have assumed the $R_+$ symmetry operator to be diagonal. A more general situation where a nondiagonal $R_+$ appears is when we do not have full knowledge of the trajectory, but only of some coarse-grained information $\mu$ at each time step (for instance, in the example of Sect. 5.1, outcome $\mu$ could be whether a $\bullet \to \bullet$, $\circ \to \circ$, $\bullet \to \circ$, or $\circ \to \bullet$ transition occurred, without knowing exactly the initial and final states of the transition). Accordingly, we write

$$P = \sum_{\mu} \Pi_{\mu}, \quad (\Pi_{\mu})_{ij} = \operatorname{Prob}(i \to j | \mu), \qquad (63)$$
where $P(i \rightarrow j|\mu)$ is the transition probability for a process $a \rightarrow b$ conditioned on the outcome $\mu$. The probability of a sequence of coarse-grained outcomes $\mu = (\mu_1, \mu_2, \dots, \mu_T)$ is given by summing over all trajectories (61) consistent with the outcomes $\mu$:

$$ \mathbb{P}(i_0 \rightarrow i_T | \boldsymbol{\mu}) = \sum_{i_1, \dots, i_{T-1}} (\Pi_{\mu_1})_{i_0 i_1} (\Pi_{\mu_2})_{i_1 i_2} \cdots (\Pi_{\mu_T})_{i_{T-1} i_T} = (\Pi_{\mu_1} \cdots \Pi_{\mu_T})_{i_0 i_T} \quad (64) $$

(Here, we assume that the exact initial and final states $i_0$ and $i_T$ are known; if that is not the case, one performs an additional sum over them.) The probability of the reverse sequence of outcomes $\mu_{\leftarrow} = (\mu_T, \mu_{T-1}, \dots, \mu_1)$ occurring is, accordingly,

$$ \mathbb{P}(i_T \rightarrow i_0 | \boldsymbol{\mu}_{\leftarrow}) = \sum_{i_1, \dots, i_{T-1}} (\Pi_{\mu_T})_{i_T i_{T-1}} \cdots (\Pi_{\mu_1})_{i_1 i_0} = (\Pi_{\mu_1}^T \cdots \Pi_{\mu_T})_{i_0 i_T} \quad (65) $$

To proceed, we assume that all $\Pi_\mu$ have the same $R_+$ symmetry, $R_+ \Pi_\mu^T R_+^{-1} = \Pi_\mu$ (naturally, $R_+$ is also a symmetry of their sum, the stochastic matrix $P$). Then, we can write

$$ \begin{aligned} \mathbb{P}(i_T \rightarrow i_0 | \boldsymbol{\mu}_{\leftarrow}) &= (R_+^{-1} \Pi_{\mu_1} R_+ R_+^{-1} \Pi_{\mu_2} R_+ \cdots R_+^{-1} \Pi_{\mu_T}^T R_+)_{i_T i_0} \\ &= \sum_{j_0, j_T} (R_+^{-1})_{i_0 j_0} (R_+)_{j_T i_T} \mathbb{P}(j_0 \rightarrow j_T | \boldsymbol{\mu}), \end{aligned} \quad (66) $$

i.e., the probability of a backward sequence of outcomes is given by the sum of the probabilities of all forward sequences compatible with $\mu$ and weighted by the entries of the $R_+$ symmetry. In the case of a diagonal $R_+$, Eq. (60), we recover a conditioned version of Eq. (62), namely

$$ \mathbb{P}(i_T \rightarrow i_0 | \boldsymbol{\mu}_{\leftarrow}) = \frac{\pi_{i_T}}{\pi_{i_0}} \mathbb{P}(i_0 \rightarrow i_T | \boldsymbol{\mu}). \quad (67) $$

## 5.6 Time Reversal of Correlation Functions

For the two symmetries that involve a minus sign, $R_-$ and $S$, the resulting operator $-L^T$ or $-L$ is not a proper generator of stochastic dynamics and, as such, the “time-reversed” trajectories are not properly defined. Nevertheless, as in the quantum case [9], certain “time-reversed” correlation functions are still meaningful. For definiteness, we consider an $S$ symmetry and a vector of observables $O$ and a probability vector $p_0$ that are both invariant under it, $OS = O$ and $Sp_0 = p_0$. The time-dependent expectation value $\langle O(t) \rangle = O \cdot p(t)$, where $p(t) = e^{Lt} p_0$ satisfies:

$$ \langle O(t) \rangle = O \cdot e^{Lt} p_0 = e^{-2\alpha t} O S \cdot e^{-Lt} S p_0 = e^{-2\alpha t} \langle O(-t) \rangle. \quad (68) $$

It follows that, even though the $S$-transformed matrix $SLS^{-1}$ does not generate Markov dynamics, the weighted expectation value $\langle \hat{O}(t) \rangle = e^{\alpha t} \langle O(t) \rangle$ is time-reversal symmetric, $\langle \hat{O}(t) \rangle = \langle \hat{O}(-t) \rangle$.

## 6 Conclusions and Outlook

In this paper, we showed that a powerful classification scheme based on (anti) unitary symmetries in quantum mechanics can be extended to classical stochastic processes. By relaxing the assumption of antiunitarity, we put forward a fifteenfold classification based entirely on linear involutions of the Markov generator $L$ of the stochastic process. However, the existence of a preferred basis in which the off-diagonal matrix elements of $L$ are positive, strongly
constrains the construction of analytical examples, a task that is usually straightforward in the quantum case. Because we were only able to parametrize L analytically in five of the fifteen classes, we then turned to a numerical optimization algorithm to find examples in the remaining classes. To this end, we minimized the sum of negative off-diagonal elements under the constraint of the symmetry transformations, finding examples in six additional classes (either reaching a zero cost function in a finite number of steps or decreasing it asymptotically). Therefore, in total, we found representative examples in eleven of the fifteen classes. While we were not able to find any examples in the last four classes, neither do our approaches—both analytical or numerical—disprove their existence. As such, we leave proving no-go theorems for these classes or explicitly building examples in them as an open problem that deserves further investigation. We also defer the improvement of the Monte Carlo optimization algorithm such that system sizes $N \gg 8$ can be reached efficiently to the future.

Our symmetry-based considerations allowed us to identify classical stochastic processes with so far unidentified properties, usually attributed solely to quantum processes, namely, the existence of dihedral spectral symmetry or Kramers' degeneracy. Furthermore, the former leads to a time-reversal-like property of expectation values of classical observables, even if the dynamics are out of equilibrium. Finally, we also interpreted $R_+$ as a generalized detailed balance and discussed its consequences for the stationary measure and stochastic trajectories of the associated Markov process. Our work opens the door for identifying these, and other unexplored, properties in realistic physical contexts.

## Appendix A Steady-State Blocks of Reducible Markov Matrices

In this appendix, we give a proof, with only a mild assumption, of the statement made in the main text that, for a reducible Markov matrix with a commuting symmetry, the block containing the steady state can also be reduced to Markov form, i.e., satisfies conditions (2)–(4). We give evidence for the accuracy of the assumption behind the proof by working out in detail an example of a stochastic process with SU(2) symmetry, where it is naturally satisfied.

We consider a reducible Markov matrix with states $\{|i\rangle\}_{i=1}^N$ and a commuting symmetry $LU = UL$. Let $|p_\infty\rangle$ be the (generically unique) steady state, $\langle 1| = \sum_i \langle i|$ the associated left eigenvector (with $\langle 1|p_\infty\rangle = 1$), and $Q_u$ be the projector into the steady-state block where $U$ has eigenvalue $u$, $Q_u U Q_u = u Q_u$; we have $Q_u |p_\infty\rangle = |p_\infty\rangle$ and $\langle 1|Q_u = \langle 1|$. Now, let $\{|u_a\rangle\}_{a=1}^d$ be a complete basis of the steady-state block, where $d$ is its dimension, i.e., $Q_u = \sum_a |u_a\rangle\langle u_a|$, and define a new basis

$$|a\rangle = \frac{|u_a\rangle}{\langle 1|u_a\rangle} \quad \text{and} \quad \langle \tilde{a}| = \langle 1|u_a\rangle\langle u_a|, \qquad (A1)$$

which also satisfies $\sum_a \langle \tilde{a} | = \langle 1 |$ and $Q_u = \sum_a |a\rangle \langle \tilde{a} |$. We assume a basis $|u_a\rangle$ exists such that $\langle i | u_a \rangle$ are real and positive. In this basis, the projection $L_u = Q_u L Q_u$ into the steady-state block is given by

$$(L_u)_{ab} = \langle \tilde{a} | L | b \rangle = \sum_{ij} \langle \tilde{a} | i \rangle L_{ij} \langle j | b \rangle = \sum_{i,j \neq i} (p_{ai} L_{ij} p_{j|i} - p_{ai} L_{ji} p_{i|i}), \quad (A2)$$

where we used the fact that for the full Markov matrix in the original basis $L_{ij} = -\sum_{j \neq i} L_{ji}$ and defined the conditional probabilities $p_{ai} := \langle \tilde{a} | i \rangle$ and $p_{i|a} := \langle i | a \rangle$. Since $\langle 1 | = \sum_i \langle i |$ and $\langle 1 | i \rangle = 1$, we have $\sum_a p_{ai} = 1$ and $\sum_i p_{i|a} = 1$. From this, it follows that $p_{ai} \in [0, 1]$.
At this point, $(L_u)_{ab}$ is not guaranteed to be of Markov form. To proceed, we assume that the steady-state block corresponds to the symmetric sector of $U$ and, furthermore, that

$$p_{a|i} \in \{0, 1\}. \qquad (A3)$$

To motivate assumption (A3), we consider three coupled two-level systems with (original) basis

$$\{|i\rangle\}_{i=1}^{8} = \{|\uparrow\uparrow\rangle, |\uparrow\downarrow\rangle, \dots, |\downarrow\downarrow\rangle\}. \qquad (A4)$$

If there is a $\mathbb{Z}_2$ spin-flip symmetry, the symmetric sector has basis

$$\{|\tilde{a}\rangle\}_{a=1}^{4} = \{|\uparrow\uparrow\rangle + |\downarrow\downarrow\rangle, |\downarrow\uparrow\rangle + |\uparrow\downarrow\rangle, |\uparrow\downarrow\rangle + |\downarrow\uparrow\rangle, |\uparrow\uparrow\rangle + |\downarrow\downarrow\rangle\}; (A5)$$

If, instead, the system has SU(2) symmetry, the symmetric sector of maximum spin is generated by

$$\{|\tilde{a}\rangle\}_{a=1}^{4} = \{|\uparrow\uparrow\rangle, |\uparrow\downarrow\rangle + |\downarrow\uparrow\rangle, |\uparrow\downarrow\rangle + |\downarrow\downarrow\rangle, |\downarrow\uparrow\rangle + |\uparrow\downarrow\rangle\}; (A6)$$

both bases satisfy assumption (A3). In the next section, we work out the SU(2) example in detail.

Equipped with assumption (A3), we can show that $L$ projected into the steady-state block, Eq. (A2), indeed satisfies the conditions (2)–(4) of a Markov matrix. First, because $p_{a|i}$, $p_{i|a}$, and $L_{ij}$ are all real, it immediately follows that $(L_u)_{ab}$ is also real. Next, to check the positivity of the nondiagonal entries of $L_u$, we note that if $a \neq b$, the last term in the right-hand side of Eq. (A2) vanishes by assumption. Since $p_{a|i}$, $p_{i|a}$ and $L_{ij}$ with $i \neq j$ are all nonnegative, it follows that

$$(L_u)_{ab} = \sum_{i, j \neq i} p_{a|i} L_{ij} p_{j|b} \geq 0, \quad \text{if } a \neq b. \qquad (A7)$$

Finally, using $\sum_a p_{a|i} = 1$, we have

$$\sum_a (L_u)_{ab} = \sum_{i, j \neq i} (L_{ij} p_{j|b} - L_{ji} p_{i|b}) = 0, \qquad (A8)$$

thus showing that probability is conserved, which concludes the proof of Eqs. (2)–(4) for $L_u$.

## A.1 Example: SU(2) Symmetry

Consider the following Markov process for the evolution of the probability vector $|p\rangle = \sum_{i_1 \dots i_N} p_{i_1, \dots, i_N} |i_1, \dots, i_N\rangle$ over the states $|i_1 \dots i_N\rangle$, $i_n \in \{-1/2, 1/2\}$:

$$L = \gamma_+ \sum_i \left[ \sigma_i^+ + \frac{1}{2} (\sigma_i^z - 1) \right] + \gamma_- \sum_i \left[ \sigma_i^- - \frac{1}{2} (\sigma_i^z + 1) \right] \\ + \gamma_0 \sum_{ij} \left[ \sigma_i^+ \sigma_j^- - \frac{1}{4} \left( \sigma_i^z + \sigma_j^z + \sigma_i^z \sigma_j^z - \frac{N(N+2)}{N^2} \right) \right], \qquad (A9)$$

where $\sigma_i^\alpha$, $\alpha = x, y, z$ are the usual Pauli matrices supported on site $j$, which describes classical spin injection and removal at each site at rates $\gamma_+$ and $\gamma_-$, respectively, and spin exchange between sites $i$ and $j$ at rate $\gamma_0$. Introducing collective spin variables $S^\alpha = \sum_i \sigma_i^\alpha / 2$, $L$ can
be written as

$$L = \gamma_+ \left( S^+ + S^z - \frac{N}{2} \right) + \gamma_- \left( S^- - S^z - \frac{N}{2} \right) \\ + \gamma_0 \left[ S^+ S^- - S^z - (S^z)^2 + \frac{N}{2} \left( \frac{N}{2} + 1 \right) \right]. \qquad (A10)$$

$L$ commutes with $S^2 = (S^x)^2 + (S^y)^2 + (S^z)^2$ and, thus, can be decomposed into $\text{su}(2)$ irreducible sectors. The steady state belongs to the symmetric sector of maximal spin, where $S^2 = s(s+1)$ with $s = N/2$, which is spanned by the $2s+1$ states ($m \in \{-2s, \dots, 2s\}$)

$$|s, m\rangle \propto \sum_{\sigma \in S_{2s}} |i_{\sigma(1)}, \dots, i_{\sigma(N)}\rangle_m, \qquad (A11)$$

with $|i_1, \dots, i_N\rangle_m = \underbrace{|1/2, 1/2, 1/2, \dots, \dots, \dots|}_{s+m} \underbrace{-1/2, -1/2, -1/2, \dots, \dots, \dots}_{s-m}$, and where the sum has

$(2s \quad s+m)$ terms, i.e.,

$$\langle \mathbf{1} | s, m \rangle = \sqrt{\binom{2s}{s+m}}. \qquad (A12)$$

Therefore, we can define

$$\langle \tilde{m} | = \sqrt{\binom{2s}{s+m}} \langle s, m | \quad \text{and} \quad |m\rangle = \frac{1}{\sqrt{\binom{2s}{s+m}}} |s, m\rangle. \qquad (A13)$$

Using this basis, we find that the Markov matrix reduced to the symmetric subspace, $S^2 = s(s+1)$, is given by

$$\begin{aligned} [L_{s(s+1)}]_{m,m'} &= \langle \tilde{m} | L | m' \rangle \\ &= \gamma_+ \left[ (s-m) \delta_{m,m'+1} - (s-m) \delta_{m,m'} \right] \\ &\quad + \gamma_- \left[ (s+m) \delta_{m,m'-1} - (s+m) \delta_{m,m'} \right], \end{aligned} \qquad (A14)$$

which is clearly a Markov matrix.

Let us illustrate these results for the specific case $N=3$. The $\langle \tilde{m} |$ basis for the symmetric sector is

$$\begin{aligned} \langle \tilde{m} &= 3/2 | = \langle 1/2, 1/2, 1/2 | , \\ \langle \tilde{m} &= 1/2 | = \langle -1/2, 1/2, 1/2 | + \langle 1/2, -1/2, 1/2 | + \langle 1/2, 1/2, -1/2 | , \\ \langle \tilde{m} &= -1/2 | = \langle 1/2, -1/2, -1/2 | + \langle -1/2, 1/2, -1/2 | + \langle -1/2, -1/2, 1/2 | , \\ \langle \tilde{m} &= -3/2 | = \langle -1/2, -1/2, -1/2 | . \end{aligned} \qquad (A15)$$

and $|m\rangle = |\tilde{m}\rangle / \langle \mathbf{1} | \tilde{m}$. An example of the conditional probabilities is

$$\begin{aligned} P_{i} &= \binom{1}{2}, \binom{-1}{2}, \binom{1}{2} | m = \frac{1}{2} = \left\langle \frac{1}{2}, -\frac{1}{2}, \frac{1}{2} \middle| \frac{1}{2} \right\rangle = \frac{1}{3}, \\ P_{m} &= \binom{1}{2} | i = \binom{1}{2}, -\binom{-1}{2}, \binom{1}{2} \rangle = \left\langle \frac{1}{2} \middle| \frac{1}{2}, -\frac{1}{2}, \frac{1}{2} \right\rangle = 1. \end{aligned} \qquad (A16)$$
Finally, the irreducible Markov matrix restricted to the steady-state sector is, thus,

$$L_{s(s+1)} = \begin{pmatrix} -3\gamma_{+} & \gamma_{-} & 0 & 0 \\ 3\gamma_{+} & -\gamma_{-} - 2\gamma_{+} & 2\gamma_{-} & 0 \\ 0 & 2\gamma_{+} & -2\gamma_{-} - \gamma_{+} & 3\gamma_{-} \\ 0 & 0 & \gamma_{+} & -3\gamma_{-} \end{pmatrix}. \quad (A17)$$

## Appendix B Explicit Examples of Markov Matrices in Different Classes

In this appendix, we write down one specific example of a Markov matrix belonging to each of the classes obtained numerically with the stochastic optimization algorithm in Sect. 4 for which no analytical solution exists. Recall that for classes $AI_-$, $CI$, and $CI_-$, the cost function reaches 0 in a finite number of steps and, consequently, the matrices displayed below for those classes satisfy the Markov property to machine precision; for classes $DIII^†$, $CI_+$, and $CI_-$, $f \to 0$ only asymptotically and, therefore, there are some residual nondiagonal negative entries in the matrices obtained. (For conciseness, all entries have been rounded to four decimal places.)

### B.1 Class $AI_-$

$$L = \begin{pmatrix} -0.2911 & 0.0004 & 0.0011 & 0.0061 & 0.1413 & 0.0038 & 0.0882 & 0.0748 \\ 0.0167 & -0.2849 & 0.0587 & 0.0577 & 0.0246 & 0.0973 & 0.1081 & 0.0902 \\ 0.0538 & 0.0058 & -0.2157 & 0.0749 & 0.0004 & 0.0457 & 0.0040 & 0.0359 \\ 0.0118 & 0.0471 & 0.0655 & -0.1896 & 0.0419 & 0.0479 & 0.0130 & 0.1512 \\ 0.0382 & 0.0011 & 0.0012 & 0.0025 & -0.2928 & 0.0058 & 0.0026 & 0.0023 \\ 0.0010 & 0.0264 & 0.0664 & 0.0001 & 0.0838 & -0.3096 & 0.0976 & 0.1434 \\ 0.0008 & 0.0716 & 0.0134 & 0.0030 & 0.0006 & 0.0716 & -0.3142 & 0.0039 \\ 0.1688 & 0.1324 & 0.0092 & 0.0453 & 0.0001 & 0.0374 & 0.0006 & -0.5016 \end{pmatrix} \quad (B18)$$

### B.2 Class $DIII^†$

$$L = \begin{pmatrix} -0.2878 & 0.0004 & 0.0379 & 0.0048 & -0.0000 & 0.0381 & 0.0583 & 0.0292 \\ 0.0000 & -0.3474 & 0.0000 & 0.0000 & -0.0000 & 0.0000 & 0.0057 & 0.0112 \\ 0.0320 & 0.0344 & -0.2938 & 0.0017 & -0.0005 & 0.0320 & 0.0161 & 0.0003 \\ 0.0000 & 0.0056 & 0.0000 & -0.3157 & 0.0140 & -0.0000 & 0.0068 & 0.0090 \\ 0.0000 & 0.1883 & 0.0000 & 0.2985 & -0.0133 & -0.0000 & 0.0395 & 0.1655 \\ 0.2557 & 0.1095 & 0.2557 & 0.0107 & -0.0002 & -0.0702 & 0.2172 & 0.1193 \\ 0.0000 & 0.0001 & 0.0002 & 0.0000 & -0.0000 & 0.0000 & -0.3528 & 0.0000 \\ 0.0000 & 0.0090 & 0.0000 & 0.0000 & -0.0000 & 0.0000 & 0.0093 & -0.3345 \end{pmatrix} \quad (B19)$$
## B.3 Class CI

$$ L = \begin{pmatrix} -0.2780 & 0.0001 & 0.0488 & 0.0008 & 0.0395 & 0.0005 & 0.1236 & 0.1520 \\ 0.0019 & -0.2660 & 0.1728 & 0.0006 & 0.0003 & 0.0003 & 0.0499 & 0.0177 \\ 0.0446 & 0.0360 & -0.4147 & 0.0010 & 0.0006 & 0.2001 & 0.0006 & 0.0001 \\ 0.0068 & 0.0229 & 0.0023 & -0.1649 & 0.0059 & 0.0016 & 0.0256 & 0.0851 \\ 0.0028 & 0.0001 & 0.0364 & 0.0009 & -0.1926 & 0.0000 & 0.0312 & 0.0056 \\ 0.0566 & 0.0866 & 0.1449 & 0.0364 & 0.0016 & -0.2665 & 0.0811 & 0.1645 \\ 0.1644 & 0.0636 & 0.0001 & 0.0006 & 0.1072 & 0.0031 & -0.3177 & 0.0009 \\ 0.0010 & 0.0569 & 0.0094 & 0.1246 & 0.0375 & 0.0608 & 0.0057 & -0.4259 \end{pmatrix} \quad (B20) $$

## B.4 Class CI++

$$ L = \begin{pmatrix} -0.2546 & 0.0546 & 0.0118 & 0.0000 & 0.0712 & 0.0154 & 0.0283 & 0.0005 \\ 0.0755 & -0.2116 & 0.0612 & -0.0000 & 0.0345 & 0.0741 & 0.0676 & 0.1150 \\ 0.0477 & 0.0013 & -0.4242 & 0.0000 & 0.0197 & 0.0770 & 0.1564 & 0.0111 \\ 0.0000 & -0.0000 & -0.0000 & -0.0000 & -0.0000 & 0.0000 & -0.0000 & 0.0000 \\ 0.0395 & 0.0099 & 0.0388 & 0.0000 & -0.3768 & 0.0428 & 0.0428 & 0.1743 \\ 0.0733 & 0.0004 & 0.1068 & 0.0000 & 0.0659 & -0.2469 & 0.0054 & 0.0493 \\ 0.0160 & 0.1382 & 0.1633 & 0.0000 & 0.0208 & 0.0006 & -0.3474 & 0.0586 \\ 0.0026 & 0.0070 & 0.0423 & 0.0000 & 0.1648 & 0.0370 & 0.0469 & -0.4090 \end{pmatrix} \quad (B21) $$

## B.5 Class CI--

$$ L = \begin{pmatrix} -0.3648 & 0.0452 & 0.0588 & 0.2540 & 0.0665 & 0.0359 & 0.0630 & 0.0580 \\ 0.0293 & -0.2365 & 0.0469 & 0.0174 & 0.0391 & 0.0531 & 0.0405 & 0.0403 \\ 0.0681 & 0.0542 & -0.1846 & 0.0791 & 0.0137 & 0.0570 & 0.0654 & 0.0248 \\ 0.1648 & 0.0002 & 0.0012 & -0.4525 & 0.0618 & 0.0507 & 0.0037 & 0.0424 \\ 0.0795 & 0.0050 & 0.0337 & 0.0091 & -0.3407 & 0.0278 & 0.0198 & 0.0201 \\ 0.0137 & 0.0383 & 0.0253 & 0.0438 & 0.0470 & -0.2993 & 0.0455 & 0.0720 \\ 0.0041 & 0.0660 & 0.0185 & 0.0222 & 0.0547 & 0.0146 & -0.2710 & 0.0551 \\ 0.0053 & 0.0276 & 0.0002 & 0.0269 & 0.0579 & 0.0602 & 0.0331 & -0.3126 \end{pmatrix} \quad (B22) $$
## B.6 Class CI<sub>-+</sub>

$$ L = \begin{pmatrix} -0.0000 & 0.0000 & -0.0000 & -0.0000 & -0.0000 & 0.0000 & 0.0000 & -0.0000 \\ -0.0000 & -0.3941 & 0.0497 & 0.0433 & 0.0516 & 0.1812 & 0.0268 & 0.0415 \\ -0.0000 & 0.0544 & -0.4025 & 0.0222 & 0.0544 & 0.0581 & 0.1834 & 0.0300 \\ 0.0000 & 0.0348 & 0.0151 & -0.2500 & 0.0561 & 0.0423 & 0.0494 & 0.0523 \\ 0.0000 & 0.0659 & 0.0533 & 0.0322 & -0.2447 & 0.0302 & 0.0242 & 0.0390 \\ 0.0000 & 0.1789 & 0.0540 & 0.0535 & 0.0168 & -0.3941 & 0.0367 & 0.0542 \\ 0.0000 & 0.0296 & 0.1835 & 0.0579 & 0.0309 & 0.0381 & -0.3726 & 0.0327 \\ -0.0000 & 0.0305 & 0.0469 & 0.0410 & 0.0350 & 0.0442 & 0.0521 & -0.2496 \end{pmatrix} \quad (B23) $$

**Acknowledgements** L.S. was supported by a Research Fellowship from the Royal Commission for the Exhibition of 1851. This work was supported by Fundação para a Ciência e a Tecnologia (FCT-Portugal) through Grant No. UID/CTM/04540/2019 (P.R.) and by the QuantERA II Project DQUANT funded through the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 101017733. T.P. acknowledges support from the European Research Council (ERC) through Advanced grant QUEST (Grant Agreement No. 101096208), and the Slovenian Research and Innovation Agency (ARIS) through the Program P1-0402. D.B.'s contribution was in part supported by the CNRS, the ENS, and the ANR project ESQuisses under contract number ANR-20-CE47-0014-01.

**Data Availability** Data generated during the current study is available from the corresponding author upon reasonable request.

## Declarations

**Conflict of interest** The authors have no conflict of interest to declare that are relevant to the content of this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Akemann, G., Baik, J., Di Francesco, P.: The Oxford Handbook of Random Matrix Theory. Oxford University Press, Oxford (2011)

2. Verbaarschot, J.J.M., Zahed, I.: Spectral density of the QCD Dirac operator near zero virtuality. Phys. Rev. Lett. 70, 3852–3855 (1993). https://doi.org/10.1103/PhysRevLett.70.3852

3. Verbaarschot, J.: Spectrum of the QCD Dirac operator and chiral random matrix theory. Phys. Rev. Lett. 72, 2531–2533 (1994). https://doi.org/10.1103/PhysRevLett.72.2531

4. Altland, A., Zirnbauer, M.R.: Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures. Phys. Rev. B 55, 1142–1161 (1997). https://doi.org/10.1103/PhysRevB.55.1142

5. Bernard, D., LeClair, A.: A classification of non-hermitian random matrices. In: Cappelli, A., Mussardo, G. (eds.) Statistical Field Theories. NATO Science Series, vol 73, pp. 207–214. Springer, Dordrecht (2002). https://doi.org/10.1007/978-94-010-0514-2_19

6. Kawabata, K., Shiozaki, K., Ueda, M., Sato, M.: Symmetry and topology in non-Hermitian physics. Phys. Rev. X 9, 041015 (2019). https://doi.org/10.1103/PhysRevX.9.041015
7. García-García, A.M., Sá, L., Verbaarschot, J.J.M.: Symmetry classification and universality in non-Hermitian many-body quantum chaos by the Sachdev-Ye-Kitaev model. Phys. Rev. X **12**, 021040 (2022). https://doi.org/10.1103/PhysRevX.12.021040

8. García-García, A.M., Sá, L., Verbaarschot, J.J.M., Yin, C.: Toward a classification of PT-symmetric quantum systems: from dissipative dynamics to topology and wormholes. Phys. Rev. D **109**, 105017 (2024). https://doi.org/10.1103/PhysRevD.109.105017

9. Sá, L., Ribeiro, P., Prosen, T.: Symmetry classification of many-body Lindbladians: tenfold way and beyond. Phys. Rev. X **13**, 031019 (2023). https://doi.org/10.1103/PhysRevX.13.031019

10. Kawabata, K., Kulkarni, A., Li, J., Numasawa, T., Ryu, S.: Symmetry of open quantum systems: classification of dissipative quantum chaos. PRX Quantum **4**, 030328 (2023). https://doi.org/10.1103/PRXQuantum.4.030328

11. Altland, A., Fleischhauer, M., Diehl, S.: Symmetry classes of open fermionic quantum matter. Phys. Rev. X **11**, 021037 (2021). https://doi.org/10.1103/PhysRevX.11.021037

12. Lieu, S., McGinley, M., Cooper, N.R.: Tenfold way for quadratic Lindbladians. Phys. Rev. Lett. **124**, 040401 (2020). https://doi.org/10.1103/PhysRevLett.124.040401

13. Kawasaki, M., Mochizuki, K., Obuse, H.: Topological phases protected by shifted sublattice symmetry in dissipative quantum systems. Phys. Rev. B **106**, 035408 (2022). https://doi.org/10.1103/PhysRevB.106.035408

14. Haake, F.: Quantum Signatures of Chaos. Springer, Cham (2013)

15. Magnea, U.: Random matrices beyond the Cartan classification. J. Phys. A **41**(4), 045203 (2008). https://doi.org/10.1088/1751-8113/41/4/045203

16. Prosen, T.: *PT*-Symmetric quantum Liouvillean dynamics. Phys. Rev. Lett. **109**, 090404 (2012). https://doi.org/10.1103/PhysRevLett.109.090404

17. Prosen, T.: Generic examples of *PT*-symmetric qubit (spin-1/2) Liouvillian dynamics. Phys. Rev. A **86**, 044103 (2012). https://doi.org/10.1103/PhysRevA.86.044103

18. Schütz, G., Sandow, S.: Non-Abelian symmetries of stochastic processes: derivation of correlation functions for random-vertex models and disordered-interacting-particle systems. Phys. Rev. E **49**, 2726–2741 (1994). https://doi.org/10.1103/PhysRevE.49.2726

19. Schütz, G.M.: Duality relations for asymmetric exclusion processes. J. Stat. Phys. **86**, 1265–1287 (1997). https://doi.org/10.1007/BF02183623

20. Giardina, C., Kurchan, J., Redig, F., Vafayi, K.: Duality and hidden symmetries in interacting particle systems. J. Stat. Phys. **135**, 25–55 (2009). https://doi.org/10.1007/s10955-009-9716-2

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Authors and Affiliations

Lucas Sá<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> • Pedro Ribeiro<fnref n="2" />,<fnref n="3" />,⁵ • Tomaž Prosen⁴,⁵ • Denis Bernard⁶,⁵

<figure><img src="image_4.png" /></figure>

Lucas Sá
lucas.seara.sa@tecnico.ulisboa.pt

Pedro Ribeiro
ribeiro.pedro@tecnico.ulisboa.pt

Tomaž Prosen
tomaz.prosen@fmf.uni-lj.si

Denis Bernard
denis.bernard@ens.fr

<footnote n="1">1 TCM Group, Cavendish Laboratory, University of Cambridge, JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="2">2 CeFEMA, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal</footnote>

<footnote n="3">3 Beijing Computational Science Research Center, Beijing 100193, China</footnote>
⁴ Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana, 1000 Ljubljana, Slovenia

⁵ Institute of Mathematics, Physics and Mechanics, 1000 Ljubljana, Slovenia

⁶ Laboratoire de Physique de l'École Normale Supérieure, CNRS, ENS & PSL University, Sorbonne Université, Université Paris Cité, 75005 Paris, France