Theory of classical metastability in open quantum systems

We present a general theory of classical metastability in open quantum systems. Metastability is a consequence of a large separation in timescales in the dynamics, leading to the existence of a regime when states of the system appear stationary, before eventual relaxation toward a true stationary state at much larger times. In this work, we focus on the emergence of classical metastability, i.e., when metastable states of an open quantum system with separation of timescales can be approximated as probabilistic mixtures of a finite number of states. We find that a number of classical features follow from this approximation, for the manifold of metastable states, long-time dynamics between them, and symmetries of the dynamics. Namely, those states are approximately disjoint and thus play the role of metastable phases, the relaxation toward the stationary state is approximated by a classical stochastic dynamics between them, and weak symmetries correspond to their permutations. Importantly, the classical dynamics is observed not only on average, but also at the level of individual quantum trajectories: We show that time coarse-grained continuous measurement records can be viewed as noisy classical trajectories, while their statistics can be approximated by that of the classical dynamics. Among others, this explains how first-order dynamical phase transitions arise from metastability. Finally, to verify the presence of classical metastability in a given open quantum system, we develop an efficient numerical approach that delivers the set of metastable phases together with the effective classical dynamics. Since the proximity to a first-order dissipative phase transition manifests as metastability, the theory and tools introduced in this work can be used to investigate such transitions through the metastable behavior of many-body systems of moderate sizes accessible to numerics.

We present a general theory of classical metastability in open quantum systems. Metastability is a consequence of a large separation in timescales in the dynamics, leading to the existence of a regime when states of the system appear stationary, before eventual relaxation towards a true stationary state at much larger times. In this work, we focus on the emergence of classical metastability, i.e., when metastable states of an open quantum system with separation of timescales can be approximated as probabilistic mixtures of a finite number of states. We find that a number of classical features follow from this approximation, for both the manifold of metastable states and long-time dynamics between them. Namely, those states are approximately disjoint and thus play the role of metastable phases, and the relaxation towards the stationary state is approximated by a classical stochastic dynamics between them. Importantly, the classical dynamics is observed not only on average, but also at the level of individual quantum trajectories: we show that time coarse-grained continuous measurement records can be viewed as noisy classical trajectories, while their statistics can be approximated by that of the classical dynamics. Among others, this explains how first-order dynamical phase transitions arise from metastability. Finally, in order to verify the presence of classical metastability in a given open quantum system, we develop an efficient numerical approach that delivers the set of metastable phases together with the effective classical dynamics. Since the proximity to a first-order dissipative phase transition manifests as metastability, the theory and tools introduced in this work can be used to investigate such transitions -which occur in the large size limit -through the metastable behavior of many-body systems of moderate sizes accessible to numerics. With continuing advances in the control of experimental platforms used as quantum simulators, such as ultracold atomic gases, Rydberg atoms and circuit quantumelectrodynamics [1][2][3][4][5][6][7], a broad range of nonequilibrium phenomena of open many-body quantum systems has been observed recently. Theoretical studies have progressed via the combination of methods from atomic physics, quantum optics and condensed matter, giving rise to a range of techniques including quantum jump Monte Carlo (QJMC) [8], simulations via tensor network [9], and field theoretical approaches [10][11][12].

CONTENTS
Often the focus of studies on nonequilibrium open many-body quantum systems is a phase diagram of the stationary state, and the related question of the structure of dissipative phase transitions occurring in the thermodynamic limit of infinite system size. This includes whether such systems can exhibit bistability (or multistability) of the stationary state, a topic covered both theoretically [13][14][15][16][17] and experimentally [18,19], and which order parameters are relevant for distinguishing the coexisting phases. Mean-field results often suggest multiple stationary states in the thermodynamic limit [12,20], however, more sophisticated (albeit still approximate) techniques such as variational approaches [21][22][23], infinite tensor network simulations [9] or a field-theoretical analysis [12] can still indicate a unique stationary state.
While it is unusual to see phase transitions at finite system sizes [24,25], first-order phase transitions in stationary states manifest at large enough finite system sizes [26] through the occurrence of metastability, i.e., distinct timescales in the evolution of the system statistics: classically, in the probability distribution over configuration space [27][28][29][30][31]; quantum mechanically, in the density matrix [32,33]. The statistics of such systems at long times can be understood in terms of metastable phases which generally correspond to the phases on either side of the transition being distinct from the unique stationary state for a given set of parameters. Therefore, already at a finite system size the structure of a possible first-order dissipative phase transition can be fully determined by investigating metastable states of the system [33], which is of particular importance for many-body open quantum systems, where exact methods are often limited to numerical simulations of finite systems of modest size.
Metastability can also emerge in complex relaxation towards a unique stationary state, even without a phase transition present in the thermodynamic limit. This is the case in classical kinetically constrained models [34][35][36][37][38][39] and recent open quantum generalizations of these models [40][41][42]. Here, the study of metastability can unfold the long-time dynamics responsible for the complex relaxation to the stationary state [32], with metastable phases corresponding to dynamical rather than static phases.
For classical systems with Markovian dynamics [27][28][29][30][31] and open quantum systems [32,43] described by the Lindblad formalism [44,45], metastability necessarily requires a large separation in the spectrum of the master operator governing the system evolution. This separation leads to metastable states residing in a space of a reduced dimension given by the slow eigenmodes of the master operator, and long time dynamics taking place within that space.
As discussed above, it is vital to develop approaches to efficiently uncover the structure of metastable manifolds (MMs) and long-time dynamics. In this work, we achieve this for classical metastability in open quantum systems, as follows. We define classical metastability as the case, where metastable states can be approximated as proba- The manifold of metastable states is described by coefficients c k [Eq. (4)], k = 2, ..., m, of decomposition between slow modes (dots for random initial pure states). The long-time dynamics takes place within that manifold, with the exponential decay of the coefficients towards the stationary state (red sphere) [Eqs. (3) and (5)]. Metastability can be observed experimentally as a plateau in the dynamics of observable averages (c) or two-point correlations (d) appearing during the metastable regime, τ t τ [Eqs. (8) and (9)]. Black (solid) lines show observable dynamics, blue (dashed) lines the approximation by slow modes holding after the initial relaxation, t τ , and red (solid) lines the stationary value achieved after the final relaxation, t τ . (e) Long timescales can also be observed in continuous measurement records, e.g., as intermittence in detection of quanta emitted due to jumps occurring in the system (two types shown in blue and red; gray -without associated quanta), with regimes of jump activity having a length comparable to the long-time relaxation. See Appendix A for details on the model. bility distributions over a set of m states, where m − 1 is the number of slow eigenmodes in the dynamics. We then show that the corrections to that approximation, which can be estimated using an approach based on the exact diagonalization of the master operator, play the role of a figure of merit in emergent classical properties of the manifold of metastable states and its long-time dynamics. Namely, we show that m states can be considered as distinct metastable phases, as they are approximately disjoint and orthogonal to one another, while their basins of attraction form the set of m order parameters used to distinguish them. Furthermore, we find that the longtime dynamics of the system can be approximated as a classical stochastic dynamics between the metastable phases. This holds both in the average system dynamics and in individual quantum trajectories, as obtained via individual runs of an experimental system, or from QJMC simulations [8,46]. Classical trajectories arise via coarse graining of these in time. The classical dynamics between the metastable phases is then responsible for the occurrence of intermittence [20,33,41] or dynamical heterogeneity [40,42] in quantum trajectories. Therefore, classical metastability is a phenomenon occurring not only on average, but in dynamics of individual quantum trajectories. Ale these results are also discussed in the presence of further hierarchy of relaxation timescales. Finally, while our approach does not rely on the presence of any dynamical symmetries (cf. Ref. [47]), we also show that the set of metastable phases is approximately invariant under any present symmetries. Thus, dynamical symmetries lead to approximate cycles of metastable phases and permutation symmetries of the classical longtime dynamics, and, as such, we find that any nontrivial continuous symmetries of slow eigenmodes of the dynamics preclude classical metastability.
In order to verify the classicality of metastability present in a general open quantum system, we develop an efficient numerical technique that directly identifies the set of metastable phases and the effective structure of long-time dynamics. If a dynamical symmetry is present, the verification of classicality is further simplified. Our approach relies on the ability to diagonalize the system master operator, which is usually possible only for moderate system sizes, while metastability may become prominent only for large system sizes. To mitigate this potential issue, we show that for classical metastability accompanied by intermittence or dynamical heterogeneity in quantum trajectories, metastable phases can be extracted from quantum trajectories through the use of large-deviation methods, such as the "thermodynamics of trajectories" [20,41,[48][49][50]. Therefore, there is potential to study classical metastability using QJMC simulations, which are generally feasible at quadratically larger system sizes than exact diagonalization of the generator. This paper is organized as follows. In Sec. II, we review the results of Ref. [32]. In Sec. III, we introduce the general approach to classical metastability in open quantum systems. We then discuss the resulting classical structure of the MM in Sec. IV. The effectively classical system dynamics emerging at large times is discussed in Sec. V. We refine these general results considering symmetries of the system dynamics in Sec. VI. Finally, we introduce numerical approaches to unfold the structure of classical metastability in Sec. VII. A concrete example of the ap- Here, H is the system Hamiltonian, while the jump operators J j provide coupling of the system to the surrounding environment. If the interactions between the system and the environment are associated to emissions of an energy quanta, the action of jump operators can be detected through continuous measurements [46], e.g., counting of photons emitted by atoms coupled to the vacuum electromagnetic field [20,[40][41][42]. Eq. (1) is a general dynamics of a time-homogeneous Markovian open quantum system, which arises for systems interacting weakly with an effectively memoryless environment.
Since the master operator L acts linearly on ρ, the evolution can be understood in terms of its eigenmatrices R k and their corresponding eigenvalues λ k = λ R k + i λ I k [52]. The real parts of these eigenvalues are not greater than 0, λ R k ≤ 0, as the dynamics in Eq. (1) is (completely) positive and trace preserving; we order the eigenvalues by decreasing real part λ R k . In particular, zero eigenvalues correspond to stationary states [53,54]. In this work, we assume a generic case of a unique stationary state R 1 = ρ ss . The system state at time t can then be then decomposed as where the coefficients c k ≡ Tr[L k ρ(0)] are bounded by the eigenvalues of L k , with L k being eigenmatrices of L † normalized such that Tr(L k R l ) = δ kl (there is a freedom of choice to normalize by scaling either R k or L k ). The values of these coefficients for a given physical state are closely tied, such that the corresponding linear combination of R k results in a positive matrix. We refer to L k and R k as left and right eigenmatrices (eigenmodes), respectively. Note that the trace-preservation of the dynamics implies that L 1 = 1, and thus beyond ρ ss other right eigenmatrices do not correspond to quantum states, Tr(R k ) = Tr(L 1 R k ) = 0 for k ≥ 2. The time scale τ of final relaxation to ρ ss from Eq. (2) can be seen to depend on the gap in the spectrum, τ ≈ −1/λ R 2 .

B. Spectral theory of metastability
Metastability corresponds to a large separation in the real part of the spectrum [32], −λ R m −λ R m+1 ; see Fig. 1(a). Times t after the initial relaxation τ t, correspond to the terms beyond the m-th in the sum in (2) being negligible (usually τ ≈ −1/λ R m+1 , so that e tλ k ≈ 0 for k ≥ m + 1) and the reduced expansion ρ(t) = ρ ss + m k=2 c k e tλ k R k + ..., where ... stands for negligible corrections [cf. Eq. (2)]. When the separation in the spectrum is big enough, it is possible to consider times τ t τ when decay of the remaining terms can further be neglected, e tλ R k ≈ 1 for k ≤ m (usually τ ≈ −1/λ R m ). This is the metastable regime, during which the system state is approximately stationary, i.e., metastable, as captured by where we defined P as the projection onto the low-lying eigenmodes of the master operator, which is trace and Hermiticity preserving [55]. From Eq. (4) the manifold of metastable states is fully characterized by the bounded coefficients (c 2 , ..., c m ) and thus it is (m−1)-dimensional. The MM is also convex, as a linear transformation of the convex set of initial states [see Fig. 1 At later times t τ , only the slow modes contribute to the evolution [cf. Eq. (3)]. Therefore, the dynamics towards the stationary state takes place essentially inside the MM [see Fig. 1 ρ(t) = e tL MM P[ρ(0)] + ..., (5) and is generated by [cf. Figs. 1(c) and 1(d)] L MM ≡ PLP.
Metastability can be observed in the behavior of statistical quantities such as expectation values or autocorrelations of system observables [32,33,56]. For a system observable, e.g., spin magnetization, we have where  7) is accurately captured by the effective long-time dynamics in Eq. (6). Importantly, during the metastable regime τ t τ , the observable average is approximately stationary [cf. Eq. (3)], before the final relaxation to o ss [see Fig. 1(c)], allowing for a direct observation of the metastability. This, however, requires preparation of an initial system state different from the stationary state, ρ(0) = ρ ss , something often difficult to achieve in experimental settings. Nevertheless, for the system in the stationary state, metastability can be observed as doublestep decay in the time-autocorrelation of a system observable. This is a consequence of the first measurement perturbing the stationary state, thus causing its subsequent evolution, which at time greater than τ follows the effective dynamics [cf. Eq. where O denotes the superoperator representing the measurement of the observable O on a system state [32,33]. The autocorrelation initially decays from the observable variance in the stationary state, O 2 ss − O 2 ss , to the plateau at Tr[OPO(ρ ss )] − O 2 ss in the metastable regime, and afterwards to 0 during the final relaxation [see Fig. 1(d)].

C. Quantitative theory of metastability
We now consider errors of approximations in Eqs. (3)-(5) and in Eqs. (8) and (9). We argue that the central figure of merit in the theory of metastability is given by the corrections to the stationarity during the metastable regime, i.e., the corrections in Eq. (4), = sup τ t τ e tL − P .
Here, X = Tr( √ X † X) denotes the trace norm for an operator X, while for a superoperator it denotes the norm induced by the trace norm [57].
Indeed, the corrections to the positivity of metastable states projected on the low-lying modes are defined by the distance to the set of density matrices, = sup with ρ and ρ(0) being density matrices [57] (the equality in the second line is proven in Appendix C), and can be bounded by the corrections to the stationarity in Eq. (10), by considering the distance to ρ ≡ ρ(t) with time t within metastable regime, where the inequality is obtained by exchanging the supremum and the infimum, which can only lead to a value increase (actually the norm is convex and continuous so the value is left unchanged). In the next section, we present how these corrections translate into corrections to the structure of metastable phases when the classical metastability occurs. Furthermore, beyond the metastable regime, the corrections in Eq. (5) decay exponentially, as in the leading order they can be shown to be bounded by 2C n MM , where n is an integer such that t/n belongs to metastable regime [58]. Similarly, the corrections to observable averages and correlations in Eqs. (8) and (9) are bounded by 2C n MM O max and 2C n MM O 2 max , respectively, where the max norm is the maximal singular eigenvalue of O.
Finally, we note that the corrections to the stationarity in Eq. (10) depend on the length of the chosen metastable regime, i.e., the choice of τ and τ such that for times τ t τ the truncation in Eq. (5) holds. Nevertheless, these corrections are always bounded from below by the corrections to the positivity in Eq. (11), independently from the choice of the length of the metastable regime. We note, however, that a pronounced metastable regime is a hallmark of metastability phenomenon, as thus some of our results rely on it being much longer than the initial relaxation time: the classical hierarchy of metastable phases is discussed in Sec. IV, and the correspondence of coarse-grained quantum trajectories to classical stochastic trajectories is discussed in Sec. V.

D. Dissipative phase transitions
When metastability is a consequence of approaching a first-order dissipative phase transition, we have by definition λ R k → 0 (and thus −1/λ R k → ∞) for k ≤ m. In this case, the timescales describing the final relaxation diverge, τ , τ → ∞. Therefore, the smallest time t τ considered in Eq. (10) (i.e., the beginning of the metastable regime) can be chosen arbitrarily large, so that the corrections vanish, C + , C MM → 0.

III. CLASSICAL METASTABILITY IN OPEN QUANTUM SYSTEMS
We now introduce the notion of classical metastability, by the virtue of classical approximation of the degrees of freedom present in the metastable regime. The corresponding corrections, together with the corrections to the stationarity and the positivity, determine the quality of classical approximations for the structure of the MM and for the long time-dynamics in Secs. IV-VI.

A. Definition of classical metastability
We define classical metastability to take place when any state of the system during the metastable regime τ t τ can be approximated as a probabilistic mixture of m states, ρ(t) = m l=1 p l ρ l + ..., (13) where ρ l ≥ 0 with Tr(ρ l ) = 1, l = 1, ..., m, represent the states, while p l ≥ 0 with m l=1 p l = 1 represent the probabilities that depend only on the initial system state ρ(0). We refer to ρ l as metastable phases (although their metastability is not assumed, but it is proven to follow in Sec. IV, together with their disjointness). The definition in Eq. (13) is motivated by the structure of first-order phase transitions in classical Markovian dynamics, where m disjoint stationary probability distributions constitute stable phases of the system, and the system is asymptotically found in a probabilistic mixture of those phases, with probabilities depending on the initial system configuration (we further discuss the assumption of the number m of phases in Appendix C 1). In later Secs. IV-V B we show that classical properties of metastable phases and long-time dynamics akin to those in proximity to a firstorder transition in a classical system follow as well.
Remarkably, any metastable state in classical Markovian dynamics can be approximated by a probabilistic mixture of approximately disjoint metastable phases [27,28,30], whether metastability results from proximity to a first-order phase transition, or from constrained dynamics as in glassy systems. In open quantum dynamics, for the bimodal case m = 2, any metastable state is a probabilistic mixture of two approximately disjoint metastable phases [32,33]. For higher dimensional MMs, however, the structure in general is no longer classical [32], as not only disjoint phases, but also decoherence free subspaces [59-61] and noiseless subsystems [62, 63] can be metastable, e.g., when perturbed away from a dissipative phase transition at a finite system size [53]. Therefore, it is important to be able to verify whether a MM of an open quantum system is classical as defined in Eq. (13). In this section, we introduce such a systematic approach, which we refer to as the test of classicality. For a given set of m candidate system states, the test of classicality enables one to verify the approximation of Eq. (13) and thus the classical metastability. Furthermore, it also facilitates a check of whether a given set of m initial states evolve into such metastable phases. Based on this last observation, in Sec. VII A we introduce an efficient numerical technique delivering sets of candidate states which, with the help of the test of classicality, can be postselected into metastable phases forming classical MMs.

B. Test of classicality
We first note that the definition of classicality in Eq. (13) has an elegant geometric interpretation with the MM in the space of coefficients being approximated by a simplex (see Fig. 2). When the MM is classical, the coefficients c k of a general state ρ can be approximated from Eq. (4) as m l=1 p l c (l) k , up to C +C + , where C is the correction in trace norm in Eq. (13) andC + is given in Eq. (12) [64]. Here, c (l) k = Tr(L k ρ l ) represent the metastable phases ρ l in the coefficient space [cf. Eq. (4)], Therefore, the classicality implies that coefficients of any metastable state can be approximated as a probabilistic mixture of metastable phases coefficients. Thus, the MM is approximated by a simplex in the coefficient space with vertices given by the metastable phases. For lowdimensional MMs (m ≤ 4), this can be verified visually by projecting a randomly generated set of initial conditions on their metastable states (in order to sample the MM) and checking that they are found approximately within the chosen simplex (cf. Fig. 2). Motivated by the structure of classical MMs in the coefficient space, we now introduce the test of classicality -a way of checking whether degrees of freedom describing metastable states during the metastable regime correspond, approximately, to probability distributions.
Degrees of freedom in the MM are described by the coefficients of decomposition into the eigenmodes R k , k = 1, ..., m, so that, with c 1 = 1, their number is m − 1. Motivated by Eq. (13), here we instead consider the decomposition in the new basis given by the projections of  (4)], the MM features the stationary state at (0, 0, 0) (red sphere) and is approximated by the simplex (blue lines) of m = 4 metastable phases (green spheres at the vertices). The dots represent metastable states from randomly generated pure states found inside (blue) and outside (black) the simplex. (right) Barycentric coordinates (p1,p2,p3) (withp4 = 1 − 3 l=1p l ) obtained by the transformation C [Eq. (15)] to the physical basis of metastable phases [Eq. (14)] yield probability distributions for states inside the simplex (green), while for states outside the simplex (black) the maximal distance C cl becomes the figure of merit for classical metastability [Eq. (19)]. metastable phases in Eq. (14), which is encoded by the transformation so thatρ l = m k=1 (C) kl R k . In particular, the volume of the corresponding simplex in the coefficient space is |det C|/(m−1)! [65]. The decomposition of a metastable state in this new basis [cf. Eq. (4)] is given by the barycentric coordinatesp l = m k=1 (C −1 ) lk c k of the simplex in the coefficient space, so thatp l = Tr[P l ρ(0)] with the new dual basis Whenρ l are linearly independent for l = 1, ..., m, |det C| > 0 and C is invertible so that Eq. (17) is well defined. In this case, Tr(P kρl ) = m n=1 (C −1 ) kn (C) nl = (C −1 C) kl = δ kl and the normalization of the dual basis in Eq. (17) is fixed by the traces of metastable states in Eq. (14) being 1.
Although for the barycentric coordinates we have m l=1p l = 1, and thus m l=1P l = 1, they do not in general correspond to probability distributions. Indeed, they are not all positive whenever a metastable state lies outside the simplex in the coefficient space corresponding tõ ρ l , l = 1, ..., m (see Fig. 2). Nevertheless, when barycentric coordinates are close to probability distributions, the corrections in Eq. (13) are necessarily small, as where stands for ≤ in the leading order of the corrections (see Appendix C 2). Here, p l = (p) l with p being the closest probability distribution in L1 norm top with (p) l =p l , so that p−p 1 = p 1 −1, and the maximum distance to the simplex is (cf. Fig. 2) while ρ l is chosen as the closest state toρ l , so that ρ l −ρ l ≤ C + [cf. Eq. (11)], and C MM bounds the approximation of ρ(t) by the projection on the low-lying modes [cf. Eqs. (10) and (16)]. Similarly, the average distance can be considered (see Appendix C 2). The corrections C cl can be efficiently estimated using the dual basis, wherep min l is the minimum eigenvalue ofP l in Eq. (17). Apart from being easy to compute,C cl carries the operational meaning of an upper bound on the distance of the operatorsP l to the set of POVMs (P l = P † l , P l ≥ 0, m l=1 P l = 1). Indeed, for P l ≡ [P l +max(−p min l , 0)]/(1+ C cl /2), l = 1, ..., m, the distance of the probability distribution p l ≡ Tr(P l ρ) top l , l = 1, ..., m, is bounded byC cl (see Appendix D 1).

C. Figures of merit
From Eq. (19), we obtain a criterion for verification of whether for a given set of states, the MM can be approximated as a probabilistic mixture of the corresponding metastable states [Eq. (14)]. In particular, since the norm in Eq. (18) is bounded by 2, whenever the MM is classical. Moreover, it can be shown that C cl m(C+C + ), where C is the approximation in Eq. (13) (see Appendix C 2). Since for a classical MM we have, by definition, C 1, Eq. (21) follows for a finite m. We thus conclude that Eq. (21) is a necessary and sufficient condition for classical metastability. In particular, approaching a first-order classical dissipative phase transition with a finite m requires C cl → 0. Interestingly, the bimodal case, m = 2, is always classical with C cl =C cl = C cl = 0 [32,33], but for higher dimensional MMs Eq. (20) can uncover the presence of classical metastability, see, e.g., Ref. [51]. Importantly, this condition is independent from the presence of dynamical symmetries.
In Sec. IV, we show that the metastable phases in Eq. (14) are approximately disjoint, while the operators in Eq. (17) take the role of basins of attractions. Moreover, in Sec. V, we show how the long-time dynamics towards the stationary state corresponds approximately to classical dynamics between metastable phases. Corrections in those results depend only on C MM , C + , and C cl defined in Eqs. (10), (11), and (19). Thus, they can be viewed as a complete set of figures of merit characterizing classical metastability in open quantum systems.

IV. CLASSICAL METASTABLE PHASES
In Sec. III, we introduced the definition of classical metastability of when MMs of open quantum systems can be approximated as probabilistic mixtures of a set of states. We now show that in this case, those states are necessarily metastable and constitute distinct phases of the system, in analogy to first-order phase transitions and metastability in classical Markovian systems [30].

A. Physical representation of metastable manifold
For states ρ l in Eq. (13), the distance to their projectionsρ l in Eq. (14) is bounded by (m + 1)(C +C + ), where C is the correction in Eq. (13) (see Appendix C 2). Therefore, for m(C +C + ) 1, these states are indeed metastable. In the test of classicality, however, ρ l in Eq. (18) are chosen as closest states to the projections ρ l in Eq. (14), and thus may be different to those initially considered in Eq. (13). Nevertheless, this change may lead to the increase of the corrections in Eq. (13) at most by C MM +C + + C + (see Appendix C 2). Furthermore, these states are also metastable, as their distance to their projections is bounded by 2C + . Finally, the distance to the initially considered states is bounded by (m + 1)(C +C + ) + C + . Therefore, for a finite m, metastable phases are uniquely defined up to the considered corrections.
In contrast to the eigenmodes of the master operator, the projections of metastable phases in Eq. (14) feature normalized trace, Tr(ρ l ) = 1, l = 1, ..., m, are Hermitian and approximately positive [see Sec. II A and cf. Eq. (11)]. Moreover, when the condition in Eq. (21) is fulfilled, any metastable state is approximated well by their probabilistic mixture [cf. Eq. (18) and Fig. 2]. Thus, the projections in Eq. (14) can be considered as physical basis of the MM and approximate metastable phases.
Furthermore, the dual basis operators in Eq. (17) determine the decomposition of a metastable state into the basis in Eq. (14), and as such, when the condition in Eq. (21) is fulfilled, they represent approximate basins of attraction for metastable phases (see also Appendix D 1). Importantly, via barycentric coordinates in Eq. (27), they define order parameters distinguishing metastable phases, Tr(P kρl ) = δ kl , with system observable averages being their linear combinations [cf. Eq. (8)].
Finally, probability distributions approximating barycentric coordinates are the physical representation of m-1 degrees of freedom present in the metastable regime. This follows from the fact that they approximate linear combinations of m − 1 coefficients that determine metastable states in Eq. (4).

B. Approximate disjointness of metastable phases
The metastable phases in Eqs. (13) and (18) can be shown to be approximately disjoint, that is, they describe states restricted to distinct regions of the system space.
First, in Appendix D 2, we prove that for the states closest to Eq. (14) we have where k, l = 1, ..., m and ρ l is the closest state to the projectedρ l . We have two types of corrections: C + , Eq. (11), bounds the distance betweenρ l and ρ l , and C cl bounds the distance of any metastable state from the simplex of metastable phases in the MM [cf. Eq. (19)]. Since √ ρ l is positive and normalized in the scalar product, the bound in Eq. (22) implies that the metastable phases are approximately disjoint. For the states in Eq. (13), the bound in Eq. (22) further reduces to √ 2C cl . Bounds on the scalar product and trace distance of metastable states are given in Appendixes D 2 and D 3.
Second, we can consider dividing the system space to capture supports of metastable phases. Indeed, for the subspace H l defined as the space spanned by eigenstates ofP l in Eq. (17) with eigenvalues equal or above 1/2, we have (see Appendix D 2) where ρ l is the closest state toρ l , k, l = 1, ..., m.
The bounds in Eqs. (23)-(25) also hold well forρ l in Eq. (14) as |Tr(1 H k ρ l ) − Tr(1 H kρ l )| ≤ C + , while for the states in Eq. (13) that project onρ l , they are further reduced to C cl , 1 − C cl and C cl , respectively. The bounds in Eqs. (23)-(25) support the statement that the metastable phases reside in approximately disjoint areas of the state space. For the bimodal case of open quantum dynamics, m = 2, approximate disjointness was already argued in Refs. [32,33].
Finally, we note that the subspace H l in Eqs. (23)- (25) captures not only majority of ρ l support, but by its definition also the corresponding basin of attraction, i.e., the initial states which evolve into metastable states close toρ l , i.e., with |1 −p l | 1. Indeed, in Appendix D 2, we show that We conclude that the basins of attractions are approximately disjoint as well.

C. Classical hierarchy of metastable phases
A further separation in the low-lying spectrum of the master operator in Eq. (1), i.e., −λ R m2 −λ R m2+1 , m 2 < m leads to a second metastable regime at times Eq. (4) and Ref. [66]]. In Appendix F, we show that for a large enough separation in the low-lying spectrum, metastable states during the second metastable regime form a classical MM, i.e., are mixtures of m 2 metastable phases. These m 2 metastable phases are approximately disjoint mixtures of m metastable phases of the first MM, and their supports are approximately disjoint [cf. Eqs. (22)- (25)]. Furthermore, each metastable phase of the first MM evolves approximately into a single metastable phase in the second MM, unless the second MM is not supported on that phase (the phase belongs to the decay subspace). Therefore, also the basins of attractions of metastable phases in the second MM are approximately disjoint.
These results are a direct consequence of long-time dynamics in a classical MM being well approximated by classical stochastic dynamics, which we discuss in next, as metastable states of classical stochastic dynamics are known to be mixtures of as many metastable phases as the number of low-lying modes [28,30].

V. CLASSICAL LONG-TIME DYNAMICS
The definition of classical metastability in Eq. (13) concerns the structure of metastable states. Remarkably, we now prove that the long-time dynamics is necessarily approximately classical, with stochastic jumps occurring between disjoint metastable phases. This is a consequence of the long-time relaxation towards the stationary state effectively taking place inside the MM.
A. Classical average dynamics of system and observables

Long-time dynamics
From Eq. (5) the evolution for times t τ ≡ −1/λ R m+1 effectively takes place on the MM with the effective generator L MM defined in Eq. (6). This generator can be expressed in the basis of the metastable phases [Eqs. (14) and (17)] as where k, l = 1, .., m, and thusW Eq. (15) and see Fig. 3(a)]. The dynamics of the system state within the MM is then determined by the dynamics of the barycentric coordinates where (p) l ≡ Tr[P l ρ(0)], so that P[ρ(t)] = By definition, the long-time evolution in Eq. (6) transforms the MM onto itself, see, e.g., Fig. 1(b). This does not guarantee, however, that the simplex of m metastable phases is transformed onto itself, as the evolution may cause states inside the simplex to evolve towards states outside, and thus an initial probability distribution (positive barycentric coordinates) acquiring some negative values at later times [see the inset in Fig. 3(b)]. Therefore, the dynamics generated byW is in general not positive [cf. Fig. 3(a)]. Nevertheless, as we discuss below, when the simplex of metastable phases is a good approximation for the MM in the sense of the condition in Eq. (21),W is well approximated by a generator of stochastic classical dynamics between metastable phases.

Classical generator
Dynamics generated byW conserves the probability, as from m k=1P k = 1 we have m k=1 (W) kl = 0 (cf. Appendix B). Furthermore, it can be shown to be approximately positive, withW approximated by the closest classical stochastic generator W [see Appendix E 1 and cf. Fig. 3(a)], where the norm . From Eq. (29) the normalized distance between the generators is bounded as Note that columns of W sum to 0, and then negativity of its diagonal terms follows from the positivity of the off-diagonal terms, so that dynamics generated by W is indeed positive and probability-conserving (cf. Appendix B). For the bimodal case of m = 2, the MM is always classical with C cl = 0, and thusW = W is exactly a generator of stochastic classical dynamics [33].

Classical system dynamics
We now discuss how the dynamics generated byW is approximated by the classical dynamics generated by W.
In Appendix E 1, we show it follows from Eq. (29) that e tW − e tW Therefore, for times t W 1 1/ √ C cl the effective dynamics in the MM is well approximated by the classical dynamics, as p(t) − p(t) 1 ≤ e tW − e tW 1 p 1 e tW − e tW 1 [where p(t) = e tWp ; cf. Eq. (27)]; see Fig. 3(c). This also holds true for the corresponding density matrices (see Appendix D 3).
When the approximation in Eq. (31) holds also for times t beyond relaxation timeτ in the MM (τ ≤τ ; see Appendix D 3), i.e., there exist time t such that the stationary state ρ ss described within the MM by (p ss ) k = Tr(P k ρ ss ) is well approximated by the stationary probability p ss of the classical dynamics W [cf. Fig. 3(c)]. Indeed, whereP ss denotes the projection onp ss , and the last inequality holds for t in Eq. (32), since, by definition of the relaxation time, for t τ we have P ss − e tW 1 1 (see Appendix E 1) [68]. In particular, when 1/(τ W 1 )  5)] can be understood as dynamics between metastable phases [Eq. (27)], governed by the tracepreserving generatorW [Eq. (26)], which can be approximated by a classical stochastic generator W [Eqs. (28) and (29)], which is both trace preserving and positive; here a negative transition rate fromρ2 toρ3 (marked by red cross) is put to 0. (b) The long-time dynamics in the barycentric coordinates (cf. Fig. 2): green simplex corresponds to t τ , blue to t = τ , and red to t = τ , while the stationary state is marked by red sphere. Positive dynamics corresponds to the simplex of metastable phases mapped onto itself, which requires all metastable phases to be mapped inside the simplex at all times. Here,ρ2 initially acquires a negative probabilitỹ is of a smaller order than √ C cl , time t in Eq. (32) exists for C cl small enough. For example, this is the case in proximity to a classical first-order phase transition with degeneracy of m phases lifted in the same order, as in this case 1/(τ W 1 ) is finite. Eq. (33) also holds true for distance in trace norm of the corresponding density matrices (see Appendix D 3). Finally, as a corollary of Eq. (33) the stationary probability distribution p ss of classical dynamical generator W in Eq. (28) is unique. Thus, the classical dynamics is ergodic with the average time spent in lth metastable phase equal (p ss ) l , l = 1, ..., m.
Similarly, not only the stationary state but all eigenmodes of the long-time dynamics in the MM can be approximated by those of the classical stochastic dynamics. In particular, in Appendix E 1, we discuss approximation of the pseudoinverse ofW in Eq. (26) by the pseudoinverse of W in Eq. (28), a result which plays an important role in the approximation of quantum trajectory statistics that we discuss in Sec. V B.
Finally, we note that quality of classical approximations for the structure of the long-time dynamics depends not only on the corrections C cl within the metastability regime [Eq. (19)], but also on the final relaxation timescale τ [cf. Eqs. (32) and (33) and Appendix E 1]. This is due to the fact that the approximation in Eq. (29) captures the fastest among the low-lying modes, while the final relaxation is governed by the slowest among them.

Classical observable dynamics
We now argue how the classical dynamics within the MM at times t τ can be observed in the behaviour of expectation values or autocorrelations of system observables.
After the initial relaxation, the dynamics of the average for an observable O depends on the evolution of the distribution between the metastable phases as where (õ) l = Tr(Oρ l ), l = 1, ..., m, are the averages of the observable O in the metastable phases. The first line corresponds to Eq. (8), while the second line follows from Eq. (31) introducing additional corrections bounded in the leading order by 2t W 1 √ C cl max 1≤l≤m |(õ) l |. Similarly, the autocorrelation (35) where (Õ) kl = Tr[P k O(ρ l )] (cf. Ref. [33]). The the first line corresponds to Eq. (9), while the second line follows from Eq. (31) with the additional corrections bounded by Therefore, when the condition Eq. (32) can be fulfilled, the dynamics of averages and autocorrelations of observables is effectively classical [cf. Eq. (33)]. In particular, if the measurement of an observable O does not disrupt supports of metastable phases, i.e., (Õ) kl ≈ δ kl (õ) l , k, l = 1, ..., m, the long-time dynamics leads to the decay of the autocorrelations exactly as the decay of the autocorrelation ofõ in the classical dynamics: from the ob-servableõ variance in p ss during the metastable regime, towards 0 achieved at t τ .

Hierarchy of classical long-time dynamics
In the case of a further separation in the spectrum of the master operator in Eq. (1), which corresponds to another metastable regime [66], the condition in Eq. (32) may generally not be fulfilled. This takes place e.g., in the proximity of a first-order dissipative phase transition when a perturbation away from the transition lifts the degeneracy of m stable phases in several different orders (so that C cl is of a lower order in the perturbation than 1/τ ). Nevertheless, in Appendix F, we show that when Eq. (32) is fulfilled for the time τ 2 of relaxation towards the second MM, i.e.,τ 2 1/ √ C cl , metastable states during the second metastable regime, ρ(t) with τ 2 t τ 2 , are approximated by the projection on the low-lying modes of the classical stochastic dynamics in Eq. (28). After the second metastable regime, when system states are restricted to the smaller second MM, the system dynamics of relaxation to the stationary state is approximated by classical dynamics taking place only between m 2 metastable phases of that manifold [cf. Sec. IV and Eq. (29)]. Furthermore, when Eq. (32) is fulfilled for that classical dynamics with respect to corrections to classicality in the second manifold, the system stationary state ρ ss is well approximated by the stationary distribution of that dynamics [cf. Eq. (33)]. For the complete discussion, see Appendix F.

B. Classical characteristics of quantum trajectories
In Sec. V A, we showed that the dynamics of the average system state can be approximated with the classical dynamics generated by a classical stochastic generator. Now we argue that this relation pertains to individual experimental realizations of system evolution [46], or trajectories of the system state sampled in QJMC simulations [8] -so called quantum trajectories (see Fig. 4). First, we show that statistics of quantum trajectories can be directly related to the statistics of time spent in individual metastable phases during the effectively classical dynamics in Eq. (6). Second, we argue how coarsegraining in time returns classical trajectories between metastable phases, which for metastable phases differing in activity is the mechanism behind the phenomena of intermittence [20,48] and dynamical heterogeneity [40,41]. Finally, we explain how system metastability can manifest itself as proximity to a first-order dynamical phase transition in the ensemble of quantum trajectories [48].

Statistics of quantum trajectories
Quantum trajectories describe the system state conditioned on a continuous measurement record, e.g., counting or homodyne measurement of photons emitted by the system due to action of jump operators in Eq. (1). In particular, the statistics of the total number of jumps in a quantum trajectory (total number of detected photons) is encoded by the biased or "tilted" master operator [20,48] with ln(Tr{e tLs [ρ(0)]}) being the cumulant generating function for the number K(t) of jumps that occurred until time t for quantum trajectories initialized in ρ. The rates of the asymptotic statistics are determined then by θ(s) = lim t→∞ ln(Tr{e tLs [ρ(0)]})/t, which is simply the eigenvalue of L s with the largest real part. We denote the associated (positive) eigenmatrix as ρ ss (s) and choose the normalization Tr[ρ ss (s)] = 1, in which case ρ ss (s) is the average asymptotic state of the system in trajectories with the probability biased by the factor e −sK(t) . Moreover, we have When s → 0 we thus obtain θ(s) → 0 and −dθ(s)/ds → j Tr[J † j J j ρ ss ] ≡ µ ss so that the average activity rate is determined by the stationary state.
The nonanalyticities of θ(s) can be recognized as dynamical phase transitions [48], in analogy to nonanalyticities of the free energy in equilibrium statistical mechanics. In particular, a first-order dynamical phase transition occurs at s c for which the maximal eigenvalue of L sc is not unique, so that the derivative k(s) = −dθ(s)/ds is no longer continuous, but features a jump at s c [20,41,48,49].
Similarly, statistics for measurements of homodyne current measurement and for time-integral of system observables can be considered [69-74] (see also Appendixes E 2 b and E 2 c).

Classical tilted generator
We now present our first result regarding classicality of quantum trajectories. We argue that the tilted master operator in Eq. (36) can be approximated by a tilted classical generator encoding the statistics in stochastic trajectories of the classical dynamics in Eq. (28).
The statistics of total activity [75-77] in classical dynamics is encoded by a biased or tilted classical generator [for reviews see Refs.
[78, 79]; cf. Eqs. (28) and (36)] where (J) kl ≡ (1 − δ kl )(W) kl , k, l = 1, ..., m encodes the transition rates in the classical dynamics, while (μ in ) kl ≡ δ kl [μ l + (W) ll ] withμ l ≡ j Tr(J † j J jρl ), k, l = 1, ..., m encodes the average internal activity in metastable phases (which here is assumed Poissonian distributed; cf. Appendix B). The biased dynamics L s in Eq. (36) can be considered as the perturbation of the master operator L in Eq. (1) with (e −s − 1)J , where J (ρ) ≡ j J j ρJ † j . Therefore, for bias |s| much smaller than λ R m −λ R m+1 , that is, the separation to the fast eigenmodes, the m low-lying eigenmodes of L s in Eq. (36) are approximated by PL s P [80], which we denote in the metastable phase basis as [cf. Eq. (26)] k, l = 1, ..., m, and we havẽ where (J) kl ≡ Tr[P k J (ρ l )], k, l = 1, ..., m. In Appendix E 2 d, we show thatJ can be approximated by W +μ, where (μ) kl ≡ δ klμl is the metastable phase activity, and thus, together with Eq. (29), we obtain the approximation by the biased generator of classical dynamics [81] whereμ tot ≡ m k=1 (J) kl +μ in l is the rate of the average jump number in metastable phases (see Appendix E 2 f); the right-hand-side equals the maximal eigenvalue of W s . For corrections to Eqs. (41) and (42) Here, however, the dynamics features additional contributions from non-Poissonian fluctuations in metastable phases (see Appendix E 2 a and cf. the next section).
For dynamics of classical systems with metastability, or, more generally, for the basis of metastable phases in Eq. (14) commuting with the dual basis in Eq. (17), m 2C cl + 4C + in Eq. (40) can be further reduced to 2C cl + 2m C + (see Appendix E 2 d).
In Appendixes E 2 b and E 2 c, we show that generators of statistics of homodyne current measurement and timeintegral of system observables in the presence of classical metastability can be similarly linked to generators of statistics in classical trajectories, but with respect to time-integrals of the average value in metastable phases. (e,f ) Not only asymptotically, but already for times after the initial relaxation, t τ , the average rate and the fluctuation rate of jump number K(t) (black solid) can be approximated by the constant contributionK from before the metastable regime and the contribution from the dynamics within the MM (blue dashed), with the latter approximated by the corresponding rate for classical total activity K cl (t) (green short-dashed) [cf. Eqs. (43) and (46)]. (g) Coarse-graining of jump records (top) in time gives values close to metastable phases activity (upper center) [cf. panel (c)], up to fluctuations which decrease with grain size (here δt = 0.7τ ; cf. Sec. V B 4). In turn, they capture the average activity of the conditional system state |ψ(t) (lower center; we plot running average over δt), and the metastable phase support where |ψ(t) is found [bottom; we plot ψ(t)|P l |ψ(t) , l = 1, 3].

Classical cumulants
We now discuss how the dynamics of the first and the second cumulants of the jump number, directly accessible in experiments, are governed by the classical long-time dynamics for times longer than the initial relaxation. In particular, we argue how the asymptotic activity and fluctuation rate are approximated by the total activity and total fluctuation rates in the classical dynamics.
These results establish a further correspondence between statistics of quantum and classical trajectories for times after the initial relaxation, i.e., during and after the metastable regime. We note that, even asymptotically, these results do not directly follow from Eq. Classical dynamics of first cumulant. For times after the initial relaxation, τ t, such that t W 1 1/ √ C cl , the rate of average jump number is approximated as [cf. Eqs. (27) and (31) and see Appendix E 2 e] where S is the pseudoinverse of the master operator L in Eq. (1), Q ≡ I − P is the projection on the fast-modes of the dynamics [cf. Eq. (4)], and (μ tot ) kl ≡ δ klμ tot l . Therefore, activity is approximated by the time-integral of the total activity in classical trajectories, K cl (t) , whose statistics in encoded by W s of Eq. Here, p ss is the stationary distribution of the classical dynamics W, and thus (p ss ) l is the average fraction of time spent in the metastable phaseρ l with the total activityμ tot l . The corrections are bounded by max 1≤l≤m |μ l | p ss − p ss 1 + W √ C cl [cf. Eqs. (33) and (29)].
Classical dynamics of second cumulant. For times after the initial relaxation, τ t, such that SQ W 1 with (δσ2) kl ≡ −δ kl 2Tr[J SQJ (ρ l )], k, l = 1, ..., m, so that (σ tot ) 2 ≡μ tot +δσ2 are rates of total fluctuations in metastable phases (see Appendix E 2 f), K cl (t) is the average of K cl (t) for lth metastable phase in Eq. (43), i.e., (p) k = δ kl , andK l ≡ Tr{P l J SQ[ρ(0)]}/Tr[P l ρ(0)] is the contribution to the jump number from before the metastable regime conditioned on the metastable phase that the system evolves into.
When the approximation in Eq. (46) is valid for times after the final relaxation, the asymptotic fluctuation rate [32,82] where R denotes the pseudoinverse of the classical generator W in Eq. (28) and the second equality follows by noting that J +μ in = W +μ tot and thus (J +μ in )p ss =μ tot p ss .
Other statistics. Similarly to Eqs. (43), (46), (45) and (50), the first and the second cumulants for homodyne measurements or for time-integrals of system observables can be related to the statistics in classical dynamics with respect to observables given by the corresponding averages for metastable phases (see Appendixes E 2 b and E 2 c). Furthermore, the rates of integrals of average and autocorrelations of system measurements discussed in Eqs. (34) and (35) can be approximated analogously.

Classical dynamics of quantum trajectories
For systems exhibiting metastability in the system dynamics, individual evolutions of the system over time typically exhibit intermittence (distinct periods of jump activity isolated in time) or dynamical heterogeneity (distinct periods of jump activity isolated both in time and space) in the emission measurement record or time-integral of observables [see Fig. 1(e)]. We now explain that these features of dynamics can be understood in terms of classical dynamics between metastable phases that differ in internal (global or local) jump activity (see also Refs. [32,33]). This establishes a direct relation between classical trajectories and time-coarse-grained records of continuous measurements.
Time-coarse-grained measurement records as classical trajectories. Consider course-graining in time of a record of jump counting measurement, where the activity in time bins of length δt for n = 0, 1, 2, .... We now argue that time-coarse-grained measurement records can be interpreted as classical trajectories between metastable phases when the internal activity dominates the long time dynamics, μ 1 W 1 , and δt is chosen long enough within the metastable regime, as in this case the activity typically attains only values of the internal activities in metastable phases [see Fig. 4 The distribution of the activity k(n) with respect to the initial system state in the nth time bin, that is, with respect to the conditional system state ρ cond (t) at time t = nδt, depends only on the corresponding metastable state. Indeed, for SQ δt τ , the average activity over all trajectories originating in ρ cond (t) is approximated by [cf. Eq. (43)] wherep l (n) ≡ Tr[P l ρ cond (t)] determines the metastable state. Similarly, from Eq. (46), the variance is approximated by (assuming additionally mC MM 1; see Appendix E 2 f) In particular, when the conditional state evolves into a single metastable phase P[ρ cond (t)] =ρ l , the average activity is approximated by the activityμ in l of lth metastable phase [cf. Fig. 4(e)]. Furthermore, the constant term is absent in Eq. (53), and the variance of k(n) decays inversely with the increasing time-bin length δt [cf. Fig. 4(f)]. Therefore, for long enough metastable regime, δt can be chosen so that the fluctuations (σ tot l ) 2 /δt of k(n) between measurement records become negligible, and thus in typical measurement records, the activity takes values approximately equal the averagẽ µ in l [cf. Fig. 4(g) (center)]. For a conditional state evolving into a mixture of metastable phases, however, a constant term is present in the variance. Nevertheless, it can be shown that it arises because of a multimodal distribution of the jump number in nth time bin. Namely, whenC cl 1, the distribution of the activity can be approximated, up to corrections 2mC MM + 2C cl + C cl , as a mixture of m distributions with averages equal internal activities of metastable phase and variances inversely proportional to δt, where the distribution corresponding to lth metastable phase, that is, with the averageμ in l , is observed with the probability approximatingp l (n). This is proved in Appendix E 2 h, by postselecting trajectories in terms of probability of the final state in the time bin, ρ cond (t + δt), evolving (on average) into a metastable phaseρ l , which, formally, corresponds to performing a measurement at time t+δt with a POVM that approximatesP l in Eq. (17) [cf. Fig. 4 We conclude that in typical measurement records and for long enough δt, activity k(n), n = 0, 1, 2, ..., takes (approximately, up to fluctuations decaying inversely in δt) only m valuesμ in l , l = 1, ..., m, corresponding to the internal activities of m metastable phases. For the bimodal case m = 2, see also Ref. [33].
Dynamics of time-coarse-grained measurement records as classical long-time dynamics. We now argue that transitions in coarse-grained measurement records are captured by the generator of the effective long-time dynamics [Eqs. (26) and (28)]. In particular, the effective lifetime of the lth metastable phase in coarse-grained trajectories is approximated by τ l ≡ 1/µ l = −1/(W) ll , l = 1, ..., m.
From the discussion above, for an initial state ρ, the distribution of activity k(0) can be approximated, up to small fluctuations, by a probability distribution over metastable phase activitiesμ in l , with probabilities approximated byp l = Tr(P l ρ), l = 1, ..., m. Analogously, the distribution of the activity k(n) in a later nth time bin is approximated byp k (n) = [(e δtW ) np ] k , where t = nδt, which is further approximated by [(e δtW ) np ] k , k = 1, ..., m [cf. Eq. (31); corrections can be further reduced to nC cl by considering discrete stochastic dynamics; see Appendix E 1 f]. Therefore, the transition matrix, i.e., the probability of observing k(n) ≈μ in k conditioned on the observation of the initial activity k(0) ≈μ in l , is approximated by the classical dynamics transition matrix (e δtW ) kl [or a discrete stochastic dynamics; see Appendix E 1 f]. This relation is further corroborated by Eqs. (43) and (46), with the average and variance of the integrated activity, In summary, when the metastable phases differ in internal activity which dominates the classical rates of long-time dynamics, coarse-grained quantum trajectories inherit these macroscopic properties, leading to intermittence characterized by the timescales of the effective dynamics. Similarly, if the metastable phases differ in local activity (see Appendix E 2 a), metastability manifests itself in quantum trajectories by dynamical heterogeneity. Therefore, metastability can be observed not only on average [cf. Eqs. (34) and (35)], but also in individual realizations of jump counting experiments or individual samples of QJMC simulations. Analogous arguments hold for the measurement of homodyne current (cf. Appendix E 2 b).

Classical metastability and dynamical phase transitions
Finally, we explain how classical metastability can manifest itself as proximity to a first-order dynamical phase transition in the ensemble of quantum trajectories [48], i.e., to a first-order nonanalyticity of θ(s).
Metastable phases as eigenmodes of tilted generator. For a unique stationary state ρ ss of L in Eq. (1) and the bias |s| small enough with respect to the gap −Re(λ 2 ), the maximal eigenmode of L s in Eq. (36) is simply approximated by the stationary state ρ ss (s) = ρ ss + ... by means of non-Hermitian perturbation theory [80]. Furthermore, in this regime we can approximate both the maximal eigenvalue of L s and its derivative by To understand how nonanalyticities arise in θ(s), we consider the dominant contributions toW s . As discussed earlier, in the presence of classical metastability, for a bias much smaller than λ R m −λ R m+1 , the m low-lying eigenmodes of L s are approximated by the eigenmodes ofW s in Eq. (39). Furthermore, when metastable phases feature nontrivial internal jump dynamics, i.e., differences in their internal jump activities are significantly higher than the transition rates of W between phases, as considered in Sec. V B 4, we havẽ where h s ≡ 1 − e −s and W hs encodes the classical statistics of time-integral of the classical observable defined by the average internal activityμ in , rather than activity of trajectories (cf. Appendix B). The additional corrections in Eq. (56) in comparison to Eq. (40) are h s W 1 /2 [alternatively, we can replaceμ in in Eq. (56) byμ orμ tot , which increases the corrections at most by h s W 1 /2]. In particular, the function θ(s) is approximated by the maximal eigenvalue of W hs in Eq. (56), while ρ ss (s) is approximated by the corresponding eigenmode p ss (h s ), For the bias |s| large enough [but small with respect to the gap to fast modes of the dynamics λ R m − λ R m+1 ], the contribution from W in W hs of Eq. (56) can be further neglected. In that case, m low-lying eigenmodes of L s are approximated as metastable phasesρ l in Eq. (14) and the corresponding eigenvalues θ l (s) = (e −s − 1)μ in l + .... Furthermore, in this regime also their derivatives are approximated as k l (s) = e −sμin l + ..., l = 1, ..., m. In particular, the maximal eigenmode corresponds to a metastable phase of L s with the maximum (for s < 0) or minimum activity (for s > 0) [cf. Fig. 4 See Appendix E 2 a for corrections.
Proximity to first-order dynamical phase transitions. For metastable phases differing in activity (or observable averages or homodyne current), Eq. (58) implies a sharp change in the derivative of θ(s), i.e., −k(s), close to s = 0 [see Fig. 4(c)]. This sharp change can be interpreted as the proximity to a first-order dynamical phase transition [20,41,48,49]. An analogous argument was made for the classical Markovian dynamics in Ref. [31].
A sharp change in k(s) around s = 0, implies in turn a large second derivative of θ(s) [see Fig. 4(d)]. In particular, d 2 θ(s)/ds 2 at s = 0 determines the rate of fluctuations in jump number, which can be approximated as [cf. Eqs. (50) and (56)] where (σ in ) 2 ≡μ in + δσ2 and the additional corrections are bounded by W 1 (1 + 2 max 1≤l≤m |μ l | R 1 ). The fluctuation rate is indeed large for the stationary state being a mixture of metastable phases with different activities [cf. Figs. 4(b) and 4(d)]. This is a consequence of long-time scales of the effective classical dynamics between metastable phases which govern the intermittence in emission records [20,83], and are captured by the resolvent in the second term of Eq. (60). In contrast, when the stationary state corresponds to a single metastable phase (so that Rμ in p ss ∝ Rp ss = 0), the fluctuation rate is finite as fluctuations originate inside that metastable phase alone [up to corrections of Eq. In terms of phase-transition phenomenology, the proximity of a first-order dynamical phase transition manifests itself in a multimodal distribution of a dynamical quantity (i.e., the jump number) in trajectories for times within the metastability regime, while at longer times in the coexistence, within individual trajectories, of active and inactive regimes that can be considered as dynamical phases [cf.

VI. CLASSICAL DYNAMICAL SYMMETRIES
We now discuss how dynamical symmetries, i.e., symmetries of the Lindblad master operator in Eq. (1), are reflected in classical MMs and in classical long-time dynamics. We find that dynamical symmetries correspond to permutations of metastable phases, and corresponding discrete symmetries of the long-time dynamics. The presence of a dynamical symmetry allows us to further simplify the test of classicality introduced in Sec. III, and the numerical approaches we introduce in Sec. VII. Finally, our results pave a way for understanding the role of symmetry in dissipative phase transitions (see also Ref. [47]).
By a dynamical symmetry we refer to the generator of the system dynamics L obeying a symmetry on the master operator level, where U(ρ) ≡ U ρ U † with a unitary operator U of a symmetry (see Refs. [53, 84, 85]). This is also known as a weak symmetry. As we consider a unique stationary state, we are interested in the case when the symmetry operator U is not itself conserved by the dynamics, so that in general L † (U ) = 0 (as the number of distinct stationary states is the same as the number of linearly independent conserved quantities [85]). For example, U can describe the translation symmetry in homogeneous dissipative systems with periodic boundary conditions. From Eq. (61) it follows that L is block diagonal in the operator basis of eigenmatrices of U. Therefore, the eigenmatrices of L, R k (and L k of L † ) can be simultaneously chosen as eigenmatrices of U (and U † ), in where φ k equals a difference in arguments of U eigenvalues (mod 2π) (cf. Fig. 5 and see Appendix G 1).

A. Symmetry and general metastability
We first discuss how dynamical symmetry in Eq. (61) affects the structure of a general MM and the long-time dynamics within it.

Symmetry of metastable manifolds
As the set of all density matrices is invariant under any symmetry, its image under the dynamics featuring a dynamical symmetry is also symmetric at any time t, U[e tL (ρ)] = e tL [U(ρ)]. In particular, a unique stationary state achieved asymptotically is necessarily symmetric U(ρ ss ) = ρ ss (or, in the case of degeneracy, the manifold of stationary states is invariant). Similarly, the set of system states during the metastable regime, i.e., the set of metastable states, is invariant under the symmetry U. This can be seen from the MM being determined by the projection P on the low-lying modes in Eq. (4), which in the presence of a dynamical symmetry fulfills ]. This is a direct consequence of the modes of L being eigenmatrices of the symmetry, so that the coefficients gain a phase under the symmetry,

Symmetry of long-time dynamics
The dynamical symmetry in Eq. (61) together with the symmetry of the MM in Eq. (62) yields the symmetry of the long-time dynamics in the MM in Eq. (6) as where U MM ≡ PUP; see

B. Symmetry and classical metastability
We now explain how dynamical symmetries for classical metastability necessarily correspond to discrete symmetries, i.e., permutations of metastable phases.

Approximate symmetry of metastable phases
The symmetry U in Eq. (61) transforms the projections ρ 1 , ...,ρ m in Eq. (14) into U(ρ 1 ), ..., U(ρ 1 ), which are also projections of system states [e.g., U(ρ 1 ), ..., U(ρ m ) for states in Eq. (13)]. In the space of coefficients, the symmetry transformation is unitary, and does not change distances. Therefore, as the simplex with vertices corresponding toρ 1 , ...,ρ m approximates well the MM in the space of coefficients, so does the simplex of the transformed new vertices, and thus we expect the new vertices to be close to the those of metastable phases. Indeed, it can be shown that the set of metastable phases is approximately invariant under the symmetries of the dynamics. In Appendix G 2 a, we prove the action of the symmetry on the metastable phases k, l = 1, ..., m, can be understood as an approximate permutation of metastable phases, that is, where Π is a permutation matrix and n = 1, 2, ... are powers of the transformation [86]. Therefore, from Eq. (65) we obtain thatρ l is approximately transformed into π n (l) under symmetry applied n times, U n (ρ l ) −ρ π n (l) 7C cl , where π is the permutation corresponding to Π. Similarly for ρ l being the closest state toρ l we have U n (ρ l ) − ρ π n (l) 7C cl +2C + [cf. Eq. (11)].

No continuous symmetries
We now argue that any continuous dynamical symmetry acts trivially on the low-lying modes of the master operator when metastability is classical. A continuous weak symmetry is a symmetry for a Hermitian operator G. For a small enough φ, U φ is approximated by the identity transformation, and therefore for such values of φ we have Π = I in Eq. (65) with n = 1. On the other hand, U n φ = U nφ , and thus from Eq. (65) for any φ the symmetry U φ is approximated by I. But this is only possible when U φ = I, i.e., the symmetry leaves each metastable phase invariant, otherwise the corrections, as given by the Taylor series, could accumulate beyond 7C cl 1 (see Appendix G 2 b for a formal proof). Therefore, all slow eigenmodes of the dynamics must be invariant as well. As a corollary, we obtain that any nontrivial continuous symmetry of slow eigenmodes precludes classical metastability.

Symmetric set of metastable phases
We now show that the set of metastable phases can be chosen invariant under the action of a dynamical symmetry. For a discrete symmetry, there exist a smallest non-zero integer D such that U D P = P. We then have U D = I, and thus from Eq. (65) also Π D = I. Let be π be a permutation associated with Π. For each cycle in the permutation, we choose an element l and definẽ where d l is the length of the cycle π d l (l) = l (and thus D is divisible by d l ), while for the other elements of that cycle we definẽ and denote the transformation from the eigenmodes to this basis as C [cf. Eq. (15)]. This gives a symmetric set of metastable states, for which the distance to system states is again bounded by C + in Eq. (11). Furthermore, from Eq. (65) it can be shown that ρ l −ρ l 14C cl , l = 1, ..., m, and the corrections to classicality defined in Eq. (19) can increase at most by 14C cl (see Appendix G 2 c for the proofs). Therefore, without loss of generality, the set of metastable phases can be considered invariant under the symmetry. In Sec. VII A, we show how symmetric sets of m candidate sets can be generated efficiently.

Symmetry of classical long-time dynamics
The dynamical symmetry of the long-time dynamics in Eq. (63) in the basis of the metastable phases reads (see Appendix G 3 for the proof). This reduces the number of free parameters in the effective dynamics (cf. Appendix B).

Symmetric test of classicality
Here, we show that the low-lying eigenmodes are linear combinations of plane waves over cycles of metastable phases as a consequence of the symmetry of the long-time dynamics in Eq. (71). We then use this fact to simplify the test of classicality introduced in Sec. III B.
Structure of low-lying eigenmodes. Eigenvectors of the long-time dynamics generatorW correspond directly to the low-lying eigenmodes of the master operator L, as they determine the coefficients in the basis of Eqs. (14) and (17), where k = 1, ..., m. In particular, the left eigenvector of W corresponding to the eigenmode L k is simply the vector of kth coefficient for metastable phases [cf. Eq. (15)].
In the presence of a dynamical symmetry [Eq. (71)], W diagonalizes in the same basis as the corresponding permutation Π, whose eigenvectors are plane waves over the permutation cycles. Since those plane waves correspond to the eigenmodes of U MM [cf. Eq. (72)], with j = 0, 1, ..., d l − 1, d l being the length of the considered cycle, and e i2πj/d l the corresponding symmetry eigenvalue [cf. Eqs. (66) and (67)], the low-lying modes are their linear combinations, where C U is the transformation from the basis of the low-lying eigenmodes to the basis of Eq. (73), . Importantly, C U is block diagonal in the eigenspaces of Π, so that R k and L k are only linear combinations of the eigenmatrices in Eq. (73) that fulfill e i2πj/d l = e iφ k . In particular, when the symmetry eigenvalue e iφ k is unique among low-lying spectrum modes, R k and L k are necessarily proportional to Eqs. (73a) and (73b). For a single cycle, all symmetry eigenvalues are unique, so that C U is diagonal and determined by the coefficients of the single candidate phase, which is discussed in Ref. [47]. For more cycles, this is not the case, e.g., for symmetric eigenmodes.
Symmetric test of classicality. Equation (73) forms a valid basis for any symmetric set of m candidate states that are linearly independent. Thus, in order to verify whether candidate states indeed correspond to m metastable phases, the test of classicality is necessary, even in the case of a single cycle (see Appendix G 4).
Exploiting the structure of the eigenmodes, the test of classicality can be simplified as follows. First, only coefficients of cycle representatives are needed to construct C U in Eq. (75). Second, as C U is block diagonal, in order to find the dual basis to the plane waves [Eq. (73b)], only matrices of the size of the permutation eigenspaces need to be inverted [87]. The dual basis to metastable phases in Eq. (17) can then be found by the inverse transformation to Eq. (73b), that is, with the coefficients as in Eq. (73a). Finally, to estimate the corrections to classicality as in Eq. (20), it is enough to consider the elements of the dual basis corresponding to the chosen cycle representatives,P l [Eq. (73b) averaged over j], and multiply their contribution by d l [88].

VII. UNFOLDING CLASSICAL METASTABILITY NUMERICALLY
Finally, we now introduce two numerical methods to analyze the classical metastability in open quantum systems governed by Lindblad master equations. The first approach in Sec. VII A requires diagonalizing the master operator and its low-lying eigenmodes. The second approach in Sec. VII B instead utilizes quantum trajectories with probabilities biased according to their activity.

A. Metastable phases from master operator spectrum
Here, we introduce a general approach based on the low-lying left eigenmodes of the master operator in Eq. (1), which delivers the set metastable phases when the considered metastability is classical. We also discuss its effectiveness in the presence of a hierarchy of metastabilities or a dynamical symmetry. Finally, we consider how observable averages distinguishing metastable phases (i.e., order parameters) can be utilized.

Metastable phases construction
Our approach consists of the following steps (see also Fig. 6): 1. Diagonalize L to find the left eigenmatrices L k below the gap in the spectrum, k = 2, ..., m.

Construct candidate metastable states:
-diagonalize the (rotated) eigenmatrices L k , -choose the eigenstates associated to their extreme eigenvalues as initial states for candidate metastable states, -discard repetitions in candidate metastable states -cluster in the coefficients space. ii. If the number of candidate states < m, or the corrections to the classicality in Step 3.i. are not negligible, enlarge the set of candidates obtained from Step 3. by considering a random rotation of the basis of the left eigenmatrices in Step 2.
Step 1 in the above construction provides the lowdimensional description of the MM, and, as explained in Sec. III, allows for testing the approximation of the MM as mixtures of m candidate states. We choose Hermitian L k replacing conjugate pairs of eigenmodes where e −iϕ k is an arbitrary phase.
Step 2 relies on the result in Appendix H 1, that metastable states arising from extreme eigenstates of the dynamics eigenmodes can be used to approximate metastable phases in classical MMs, as long as only a single metastable phase is close to the extreme value of the corresponding coefficient c k (which we refer to as the case without degeneracy). We then discard any repetitions in the set of candidate states (in order to treat all coefficients on equal footing, we set the normalization c max and c min k are extreme eigenvalues of L k ). Indeed, for a given left basis we obtain 2(m−1) candidate metastable states corresponding to 2(m−1) extreme eigenvalues of the basis elements, which may provide up to m metastable phases.
In the case without degeneracy, each candidate corresponds to one of m metastable phases. In the case with degeneracy, some of extreme eigenstates may correspond Figure 6. Metastable phases construction. We sketch the approach to construct a candidate set of m metastable phases within a given MM that provides the best classical approximation in terms of corrections in Eq. (19). The construction can be refined by considering dynamical symmetries.
to mixtures of metastable phases: provided that the set of candidate states features all metastable phases, such a candidate state should be discarded in Step 3 i (this relies on the result from Appendix H 3, that the simplex of metastable phase is approximately the largest simplex inside the MM; cf. Fig. 2). However, such mixtures may cause less than m candidates to remain after clustering, or result in a set of phases which provides a poor approximation to the true MM; even without degeneracy, it is possible that some metastable phases may reside on the interior of the hypercube defined by the extreme values of the coefficients, and as such will not appear in the set of candidate states taken from extreme eigenvalues of the eigenmodes. Nevertheless, random rotations in Step 3 ii ensure that each metastable phase is eventually exposed, i.e., a basis in which that metastable phase achieves an extreme coefficient value without degeneracy is eventually considered (with the probability 1 achieved exponentially in the number of rotations; cf. Appendix H 2). Naturally, instead of considering distances between candidate states in the space of coefficients, as used in Steps 2 and 3 i, candidate states can be clustered with respect to the trace distance in the space of density matrices, while the best m candidate states can be chosen to achieve minimal corrections to classicality instead of the maximal simplex volume. These modifications, however, require working with operators on the system Hilbert space, rather than on the space of coefficients, and thus are in general more expensive numerically [for classical MMs, m ≤ dim(H)], while not necessary for MMs with nonnegligible volumes (cf. Appendixes H 3 and H 2).
Our approach will not deliver a set of metastable phases with negligible corrections to classicality whenever the MM is not classical. In particular, quantum MMs [32] featuring decoherence free subspaces [59][60][61] or noiseless subsystems [62, 63]: which, e.g., at m = 4 amount to Bloch-sphere in the coefficient space, rather than a tetrahedron (see also Ref. [89]), cannot be approximated as probabilistic mixtures of m metastable phases. Even for classical metastability emerging in many-body open quantum systems, the approach relies on the condition in Eq. (21), which may be fulfilled only at larger system sizes, when the low-lying part of the master operator spectrum to be sufficiently separated from the fast modes. In this case, our approach may not succeed for smaller system sizes with less pronounced metastability.
While degeneracies are unlikely to occur in a generic model, they typically appear in the presence a hierarchy of metastabilities or dynamical symmetries. Before discussing how to remedy these cases below, we note that the method has been recently successfully applied in Ref. [51] to the open quantum East model [40], featuring both a hierarchy of two metastable regimes and the translation symmetry.

Construction for hierarchy of metastable manifolds
In the presence of hierarchy of metastabilities with a further separation at m 2 < m in the spectrum of the master operator, the simplex of m metastable phases projected onto the coefficients (c 2 , ..., c m2 ) is approximated by a simplex with m 2 vertices corresponding to m 2 metastable phases of the second classical MM. This requires (at least) m 2 metastable phases in the first MM to evolve directly into m 2 metastable phases of the second MM. Each of other m − m 2 metastable phases of the first MM either evolves into a single phase of the second MM, or it belongs to the decay subspace in which case it in general evolves into a mixture of m 2 metastable phases of the second MM. In the former case, the metastable phases in the first MM that evolve into the same phase in the second MM are degenerate in the coefficients (c 2 , ..., c m2 ). In the latter case, they do cannot take extreme values of those coefficients, even after a rotation of the first m 2 modes (see Appendixes F and H 4). Nevertheless, rotations of all m modes allow both for the degeneracy to be lifted and for every metastable phase to take an extreme value in one of the coefficients.

Construction for metastable manifolds with symmetries
In the presence of a dynamical symmetry U [Eq. (61)], the low-lying eigenmodes of the dynamics are linear combinations of the plane waves over the cycles induced by the symmetry on m metastable phases [cf. Eq. (74)]. For an eigenmode L k , k = 2, ..., m with a symmetry eigenvalue e iφ k , U † (L k ) = e iφ k L k , and the minimal integer n k > 0 leading to e in k φ k = 1, L k is supported on cycles with the length equal n k or larger but divisible by n k . The latter case, of subcycles with length n k , leads to degeneracy of the coefficient c k for the metastable phases connected by U n k [cf. Eq. (73b) and Fig. 5(a), wherẽ ρ 1 andρ 3 are degenerate in c 2 and c 4 (n 2 = n 4 = 1)]. Nevertheless, coefficient degeneracy can be remedied and the structure of the low-lying eigenmodes can be used to actually enhance the introduced approach, as we now explain (see also Appendix H 5).
We note that when GCF(n k , n l ) < n k for all n l = n k , k, l = 2, ..., m, the eigenmode L k is supported only on cycles of length n k , as there are no longer cycles with the length divisible by n k . The symmetry U n k then acts trivially on the corresponding metastable phases and the above discussed degeneracy is absent. Analogously to the general case, any plane wave in L k can then be exposed by random rotations of the eigenmodes with the same symmetry eigenvalue e iφ k , while at least a single metastable phase from each cycle can be obtained by considering extreme eigenstates of both L R k and L I k [90]. Other metastable phases in the considered cycles can be recovered by applying the symmetry n k − 1 times, so there is no need to consider eigenmodes L l with a different symmetry eigenvalue supported on the same cycles (i.e., L l with e iφ k = e iφ k but n l = n k ).
Furthermore, metastable phases in cycles with the length corresponding to subcycles, e.g., invariant metastable phases, can also be found. For L l such that n l is divides only n k > n l for the above considered eigenmodes L k , when the degeneracy of e iφ l equals the number of already considered cycles with n k divisible by n l (i.e., the sum of the corresponding symmetry eigenvalue degeneracy for all such n k values), L l is supported on the already considered cycles. Otherwise, L l features new cycles with the length n l , which can be unfolded, as before, by rotations of all eigenmodes with the same symmetry eigenvalue as e iφ l and considering both L R l and L I l [90]. Here, equal mixtures of already considered phases connected by U n k /n l will also be found, but such candidate states will not lead to the maximal volume simplex. Again, by applying symmetry U the full (sub)cycles can be recovered, and other eigenmodes L j with n j = n l , but a different symmetry eigenvalue can be discarded. Analogous results hold for the remaining eigenmodes, but with respect to L l and e iφ l degeneracy.
In summary, the set of eigenmodes considered in Step 2 is significantly reduced, with its size equal to the number of cycles and the subcycles with other cycle's length. Furthermore, only rotations of eigenmodes with the same symmetry eigenvalue are necessary in Step 3 ii. Importantly, as the eigenstates of L R k and L I k and the corresponding metastable states are invariant under U n k , that is, generate cycles of candidate states with the length dividing n k , following the prescription above, we arrive at an invariant set of candidate states. This invariance can be maintained by clustering whole cycles of candidate states rather than individual states in Step 2. Then, without loss of generality, in Step 3 i, instead of considering subsets of all candidate states, we can choose candidate states as sets of cycles with their lengths summing to m. In that case, the volume of the simplex can be efficiently calculated as |det C U |/(m − 1)! l √ d l , where C U is the block-diagonal matrix in Eq. (75) and the product, which runs over cycle representatives, is the same for all sets of linearly independent candidates [87], while the corrections to the classicality can be efficiently calculated with the symmetric test of classicality of Sec. VI B 5.

Construction utilizing order parameters
We note that instead of considering the eigenmodes of the dynamics, we can choose a left eigenbasis formed by a set of m observables O l , l = 1, ..., m projected onto the low-lying eigenmodes, i.e., (4)], provided that those projections are linearly independent. In this case, the extreme eigenstates of P † (O k ) will give metastable states attaining extreme values in the average of the observable O k . Those metastable states will correspond to metastable phases, up to degeneracy of metastable phase averages of O k (in particular, in the presence of nontrivial dynamical symmetry of low-lying eigenmodes, it is necessary to consider observables breaking the symmetry). Among others, this can be helpful when the volume of MM in the space of coefficients is negligible. In the next section, we extend this approach by considering continuous measurements instead of system observables.

B. Metastable phases from biased quantum trajectories
In some systems, the metastability can be a collective effect emerging as the system size increases [33]. If large system sizes are required for prominent metastability, exact diagonalization may not be feasible. Therefore, we now introduce an alternative numerical approach to finding metastable phases in classical MMs using QJMC simulations [8] and biased sampling in the framework of large deviation theory (see Ref. [91] for a review). In classical stochastic dynamics biased sampling can be efficiently incorporated into the generation of trajectories, with techniques such as transition path sampling [92] and cloning [93].
Trajectories of the biased master equation L s in Eq. (36) can be viewed as trajectories of L with their probability multiplied by e −sK(t) , where K(t) is the total number of jumps K(t) occurring in a quantum trajectory of length t. The maximal eigenmode ρ ss (s) of L s corresponds then to the asymptotic system state in quantum trajectories averaged with the biased probability. In Sec. V B 5, we argued that ρ ss (s) can approximate metastable phases of extreme activity for appropriately chosen s when the activity dominates the transition rates of long-time dynamics [cf. Figs. 4(b) and 4(c)]. If the efficient biased sampling could be generalized to QJMC sampling, metastable phases with extreme activity could be accessed via time-average of a biased trajectory over time-length within the metastable regime [94].
Similarly as in the case of degeneracy of coefficients, when more than a single metastable phase corresponds to the maximum or the minimum activity, ρ ss (s) corresponds to a mixture of the metastable phases with the extreme value (e.g., when both L and J obey a symmetry that is broken in the MM, the mixture is symmetric). Nevertheless, the discussion in Sec. V B 5 is analogous for the activity of individual jumps, and thus a further distinction between metastable phases can be enabled this way, e.g., by breaking the translation symmetry of L s in the case of identical local jumps (see Appendix E 2 a).

VIII. CONCLUSIONS
In this paper, we formulated a comprehensive theory for the emergence of classical metastability for open quantum systems whose dynamics is governed by a Lindblad master operator. We showed that classical metastability is characterized by the approximation of metastable states as probabilistic mixtures of m metastable phases, where m is the number of low-lying modes of the master operator. Namely, in terms of the corresponding corrections, metastable phases are approximately disjoint, while the long-time dynamics, both on average and in individual quantum trajectories, is approximately governed by an effective classical stochastic generator. Furthermore, any nontrivial dynamical symmetries present at long times are necessarily discrete as they correspond to approximate permutations of metastable phases, under which the classical dynamics is invariant.
To investigate metastability for a given open quantum system, we introduced the test of classicality -an approach to verify the approximation of the MM by a set of candidate metastable phases. We also developed a complementary numerical approach to deliver sets of candidate metastable phases. Since that approach requires diagonalization of the master operator -a difficult task in systems of larger size -we also discussed an alternative based on the concept of biased trajectory sampling.
The introduced set of techniques allows us to achieve a complete understanding of classical metastability emerging in an open quantum system. A concrete application of the methods described here to a many-body system is found in Ref. [51], where the long-times dynamics of the glassy open quantum East model [40] is studied in detail.
An important question is how to investigate systems exhibiting general quantum metastability, e.g., featuring metastable coherences [32]. Extending the methods described here to those systems would inform the study of a general structure of first-order dissipative phase transitions in open quantum systems. We acknowledge financial support from EPSRC Grant no. EP/R04421X/1, from the H2020-FETPROACT-2014 Grant No. 640378 (RYSQ), and by University of Nottingham grant no. FiF1/3. I.L. acknowledges support by the "Wissenschaftler-Rückkehrprogramm GSO/CZS" of the Carl-Zeiss-Stiftung and the German Scholars Organization e.V., and through the European Union's H2020 research and innovation programme [Grant Agreement no. 800942(ErBeStA)]. We are grateful for access to the University of Nottingham Augusta HPC service.
[52] It is possible that the master operator may not be completely diagonalizable, in which this case its simplest form is the Jordan normal form. Under exponentiation, this gives rise to a polynomial dependence on the time, however this is accompanied by the usual exponential evolution. Our results in Secs. III, IV, V A, VI, and VII carry directly to that case by assuming that after the initial relaxation, any polynomial evolution of fast modes is dominated by the decaying exponential and can be neglected [cf. Eq. (3) Thus the set of low-lying modes is invariant under the Hermitian conjugation. [57] We restrict the discussion to the real space of Hermitian operators. In this case, the induced norm of a superoperator can be shown to be achieved for a pure state (i.e., a rank one operator). Let X be a Hermitian operator with eigenvalues x k and projections on the corresponding eigenstates denoted as ρ k . For a superoperator Y, [58] For an integer n ≥ 1 such that τ t ≡ t/n τ , that is, t belongs to the metastable regime, we have and c min k denoting the maximal and minimal eigenvalue of L k . We then choose time t within the metastable regime for which P[ρ(0)]−ρ(t) is minimal.  1 for all prime factors n of n (see Appendix G 2 a).
[87] When momenta of plane waves over candidate state cycles do not correspond to arguments of symmetry eigenvalues for the low-lying eigenmodes, C U is not invertible and det C U = 0.
[88] This follows from the spectrum of any operator unchanged under the action of the symmetry, and U † [P π n (l) ] =P π n+1 (l) . Alternatively, the minimal eigenvalue can be found from the wave-plane basis by considering averages of

APPENDIX: PROOFS
In the following appendixes, we provide detailed proofs of the results throughout the paper.
In the limit the dynamics is in the proximity of a dissipative phase transition with four independent sectors corresponding to two decoupled pure states |1 and |2 two effective twolevel atoms formed by |3 and |5 , and by |4 and |6 (see Fig. 7). The coherences between two decoupled states are not stable because of dephasing caused by J 1 and J 2 in Eq. (A3). We note that due to a perturbation of this finite-size classical phase transition by the Hamiltonian δH = ω 1 (σ x 13 + σ x 14 ) + ω 2 (σ x 25 + σ x 26 ), the long-time dynamics arises in the second order in ω 1 and ω 2 , while the structure of the phases is perturbed already in the first order [32]. Therefore, the condition in Eq. (A4) leads to the classical metastability with m = 4 in Figs. 1-4 (for parameters see below), as the dynamics caused by δH can be neglected before and during the metastable regime, t τ . Replacing J 4 in Eq. (A3) by J † 4 introduces a discrete dynamical symmetry under the simultaneous swap |3 ↔ |4 and |5 ↔ |6 (cf. Fig. 7 and see Sec. VI). This is the model depicted in Fig. 5, in the limit of Eq. (A4).

Appendix B: Classical stochastic dynamics
Here, we review properties of classical stochastic dynamics and statistics of its trajectories. We also discuss dynamical symmetries.

Stochastic trajectories
The matrix W can be considered as a Markovian generator of stochastic trajectories of system configurations. The waiting time for a transition from the kth configuration is distributed exponentially with the rate −(W ) ll , so that the average lifetime is being a diagonal matrix of activities. The (minus) first derivative of the maximal eigenvalue determines the average activity where the diagonal matrix of the system observable The rate of the average time-integral while the rate of fluctuations of time-integral Eqs. (B11) and (B12) follow from the perturbation theory for W h with respect to hm.

Dynamical symmetries
Stochastic dynamics features dynamical symmetry Π when where Π is a permutation matrix between m system configurations. Therefore, the symmetry in Eq. (B13) is equivalent to where k, l = 1, ..m and π denotes the permutation corresponding to Π. From Eq. (B14), all configurations that belong to the same cycle, feature the same decay rate, (W) kk = (W) π n (k)π n (k) , n = 1, 2, .... Furthermore, consider two different cycles of length m 1 and m 2 and let m 1,2 be the greatest common factor of m 1 and m 2 . Transitions from the first to the second cycle are the same for all elements of a subcycle of length m 1 /m 1,2 and all elements of a subcycle m 2 /m 1,2 (and analogously for transitions from the second to the first cycle). In particular, the transition rates from (or to) an invariant configuration l [π(l) = l] are the same for all the elements of a cycle, (W) kl = (W) π n (k)l [(W) lk = (W) lπ n (l) ], n = 1, 2, ....

Appendix C: Classical metastability in open quantum systems
In this appendix, we first discuss the correctness of the definition of classical metastability in Eq. (13) in terms of the number of metastable phases. We then consider the test of classicality and prove that C cl in Eq. (19) is the maximal distance of the barycentric coordinates of a metastable state from a given simplex of candidate metastable phases when measured by L 1norm. We further prove the bound on the corrections in Eq. (18). We also derive a similar bound on the average distance.Finally, we consider the optimality of the metastable phase construction in Eq. (18) in the context of corrections Eq. (13). As a corollary of derivations presented in this appendix, the second line of Eq. (11) follows.

Definition of classical metastability
In Eq. (13) we assumed the number of states to be equal to the number m of low-lying modes in the spectrum of Lindblad master operator in Eq. (1). We now justify this assumption.
A higher number of states than m Eq. (13) necessarily leads to linearly dependent matrices after the projection onto the low-lying modes as in Eq. (14). Therefore, the decomposition of P[ρ(0)] into such states in the space of coefficients in not unique [even with the additional assumption of the (approximate) positivity of decomposition]. Therefore, when the states in Eq. (13) are well approximated by their projection on the low-lying modes in Eq. (14), we conclude that there are no more than m states for the decomposition to be uniquely defined [uniqueness here is defined up to corrections in Eq. (13)]. Furthermore, when some of the states in Eq. (13) differ significantly from their projection on the low-lying modes, i.e., are not metastable, not all probabilistic mixtures of the states in Eq. (13) are metastable, and thus the probabilities, although possibly uniquely defined, do not represent degrees of freedom of the MM. Moreover, although in this case the candidate states can be replaced by the states closest to their projections in Eq. (14) [with an increase of corrections in Eq. (13) by at most C MM +C + + C + ; see Eq. (C18) below], due to linear dependency in the space of coefficients, this again will lead to a non-unique decomposition.
A lower number of states than m in Eq. (13) indicates degeneracy of the description of the MM in the space of coefficients (see e.g., Ref. [30]), leading to the effective lower dimension of the metastable state manifold than m − 1. This case will be discussed elsewhere.
Note that the choice of optimal p in Eq. (C2) is generally not unique.
We conclude that the maximal distance of barycentric coordinates to the simplex of probability distributions over all initial states of the system is given by Eq. (19). As a barycentric coordinatep l is bounded from below by the minimal eigenvalue ofP l in the dual basis, from Eq. (C3) we arrive at Eq. (20).
The average distance C cl for uniformly distributed pure initial states of the system is bounded by withP − l beingP l restricted to its negative eigenvalues and dim(H) the dimension of the system Hilbert space.
Furthermore, the average corrections in Eq. (13) are bounded as From the triangle inequality we can compare the classical approximation of metastable states in the trace norm as where in the second inequality we used Eq. (10), while the third inequality follows from ρ l − ρ l ≤ C + for ρ l being the closest state toρ l [cf. Eq. (11)] and ρ l ≤ ρ l − ρ l + ρ l ≤ 1 + C + . Therefore, using Eq. (19) we arrive at Eq. (18).

c. Optimality of test of classicality
We now consider how the construction of classical approximation in Eq. (18) compares to a given set of states in Eq. (13).
Bound on C cl in terms of corrections in Eq. (13). We now bound the corrections C cl in Eq. (19) by the corrections in Eq. (13). We have where l = 1, ..., m.
In the first inequality we used Eq. (C21) below], and in the last inequality we used the fact that P l max ≤ 1 + C cl [see Eq. (D5)] and where C are the maximal corrections in Eq. (13) over the set of metastable states. Therefore, Bound on the stationarity of states in Eq. (13). We now prove that for states in Eq. (13) we have and thus states in Eq. (13) are metastable for the corrections in Eq. (13) fulfilling mC 1. Consider ρ l in Eq. (13) as an initial state of the system, l = 1, ..., m. We denote ρ l (t) the corresponding state at time t. By definition, there exists a probability distribution p (l) k , k = 1, ...m such that for all times within the metastability regime and therefore from Eq. (C8) On the other hand, we also have [cf. Eqs. (14) and (C7)] Together with Eq. (C11), we then obtain and, thus, again using Eq. (C8) arrive at Bound on corrections for metastable phases replaced by their projections. We have where in the second inequality we used Eq. (C8). Therefore, from the triangle inequality Bound on corrections for metastable phases replaced by closest states to their projections. Let ρ l be the closest state to the projection P[ρ l ] of ρ l in Eq. (13), so that P[ρ l ] − ρ l ≤ C + [cf. Eq. (11)]. From the triangle inequality we then have [cf. Eq. (C18)] Equations (C18) and (C19) actually hold for any number of states in Eq. (13) (see the discussion on the correctness of classicality definition in Appendix C 1). We conclude that corrections in Eq. (13) can increase in the leading order by at most C MM +C + + C + when the states in Eq. (13) are replaced by closest states to their projections.
The Von Neumann trace inequality. In the proofs above, we make use of the von Neumann trace inequality for operators, which reads [97-99] where λ (X) n are singular eigenvalues of X [i.e., the eigenvalues of √ X † X] ordered decreasing in value, and dim(H) is the dimension of the system Hilbert space. In particular, we have λ (X) n ≤ X max and n λ (X) n = X , so that 3. Derivation of Eq. (11) Using the results of this appendix, let us consider a projection on the MM in Eq. (4). As P preserves the trace, from Eqs. (C1)-(C3) it follows that the closest state (density matrix) diagonal in the eigenbasis of P[ρ(0)] is determined by Eq. (C2), leading to the distance in the trace norm equal P[ρ(0)] − 1 [cf. Eq. (C3)]. On the other hand, for any density matrix ρ from the triangle inequality we have and thus the choice of ρ as diagonal in the eigenbasis of P[ρ(0)] with the eigenvalues given as in Eq. (C2) is optimal.

Appendix D: Classical metastable phases
Here, we prove the bounds in Eqs. (22), (23), and (24) of the main text. We begin by discussing the properties of the dual basis, which we then use in derivations of Eqs. (22), (23), and (24). We also prove the bound on the scalar product of the metastable phases. We finish by discussing the relations between distances in the barycentric coordinates for metastable phases and in the space of density matrices.

Properties of dual basis in classical metastable manifolds
Here, we discuss properties of minimum and maximum eigenvalues of the dual basis in Eq. (17) and find corresponding bounds on its norms in terms of C cl in Eq. (19). We also show thatC cl in Eq. (20) can be understood a distance to operators in a certain POVM.
As the norm P l max is by definition equal the maximal or minus the minimal eigenvalue ofP l , we conclude

Consider operators
We have that P l ≥ 0 and m l=1 P l = 1, so that these operators constitute a POVM, i.e., a set of operators for which p l ≡ Tr(P l ρ), l = 1, ..., m, corresponds to a probability distribution for any density matrix ρ.
We now estimate the distance of this probability distribution to the barycentric coordinatesp l ≡ Tr(P l ρ) [cf. Eq. (17)]. We have where the first inequality corresponds to the triangle inequality.

Disjointness of phases in classical metastable manifolds
a. Proof of Eq. (22) Here, we find an upper bound on the scalar product between square roots of the closest states toρ l in Eq. (14), l = 1, ..., m.
We have wherep min k andp max k are minimal and maximal eigenvalues ofP k , which are introduced to obtain products with positive operators. Ifp min k ≥ 0 we instead replace it by 0 in Eq. (D8), whilep max k ≤ 1 we instead replace it by 1. From the Cauchy-Schwarz inequality we further obtain, In the third line we used (P k −p min Note that we have |p min Eq. (D2)] and P k max ≤ 1 + C cl /2 [cf. Eq. (D5)], where C cl bounds the distance from the metastable phases simplex of any metastable state [cf. Eq. (19)]. Therefore, from Eq. (D8) we arrive at where C + is a bound from above on the distance ρ l − ρ l [cf. Eq. (11)]. In the leading order Eq. (D11) gives Eq. (22). Finally, note for the states ρ l that project ontoρ l [cf. Eqs. (13) and (14)],ρ l = P(ρ l ), in Eqs. (D12)

b. Proof of Eqs. (23)-(25)
We now prove Eqs. (23)- (25). We also prove bounds on the support of general system states that allow for showing approximate disjointness of basins of attractions.
Therefore, from Eq. (D29) we obtain and from Eq. (D30) as well as k =l where we used k =l Tr(P l ρ) = 1 − Tr(P l ρ).

c. Bound on the scalar product of metastable phases
We now prove the following bound on the scalar product of the metastable phases where k = l, k, l = 1, ..., m, ρ l is the closest state tõ ρ l in Eq. (14), and ρ l max ≤ Tr(ρ 2 l ). The same inequality holds for Tr[ρ kρl ], with ρ l being the closest state toρ l in Eq. (14), while for the states in Eq. (13) we have a reduced bound C cl /2( ρ k max + ρ l max ). Note that the scalar product of metastable phases is affected by mixedness of both phases, |Tr(ρ k ρ l )| ≤ Tr(ρ 2 k ) Tr(ρ 2 l ), which enhances their approximate orthogonality.

Trace-norm vs. L1-norm in classical metastable manifolds
Here, we discuss how the distances measured in the trace norm in the space of density matrices and by L1norm in the barycentric coordinates of the MM are related.

b. Distance between metastable phases
From Eqs. (D41) and (D42) we obtain that the distance between the projections of the metastable phases in Eq. (14) is bounded as and the lower bound also holds for any states that projected on the low-lying modes give matrices in Eq. (14) [cf. Eq. (13)]. Therefore, for the case bimodal case m = 2, we obtain that the metastable phases are disjoint with respect to the trace norm (cf. Ref. [33]).

c. Relaxation time
We now connect the relaxation time towards the stationary state ρ ss in the trace norm to the relaxation of metastable state within the classical MM. We have In the first inequality we used Pρ =  (2) and (4)] is well captured by the decay of the probabilities between metastable phases. However, from Eq. (D45), the relaxation timeτ with respect to L1-norm is generally longer than that of the relaxation time τ with respect to the trace norm τ ≤τ . (D48) In the case of a perturbation away from a classical firstorder phase transition, however, both relaxation times are of the same order in the perturbation.

Appendix E: Classical long-time dynamics
Here, we prove that the long-time dynamics is effectively classical as discussed in Sec. V. First, in Appendix E 1, we consider the dynamics of the average system state within a classical MM, and derive bounds on its approximation by classical dynamics governed by a classical stochastic generator. Second, in Appendix E 2, we prove that classical trajectories of that classical stochastic generator capture statistics of quantum trajectories on timescales longer than the initial relaxation.

Classical dynamics of average system state
We now show that the long-time dynamics of the average system state is effectively classical. We first derive the best continuous approximation of the generator of the long-time dynamics in Eq. (26) by classical stochastic dynamics between metastable phases in Eq. (28). We then prove an upper bound Eq. (29) on the distance between the two generators. We use this result to show the closeness of the effective dynamics to that generated by a classical stochastic generator in Eq. (31) and of the resulting stationary states in Eq. (33). Similarly, all timescales of the dynamics in the MM can be approximated by the classical stochastic dynamics. Here, we consider here an approximation of the pseudoinverse of the long-time dynamics generator. Finally, we also discuss the approximation of the long-time dynamics by discrete, rather than continuous, positive and trace-preserving dynamics, and prove results analogous to Eqs. (29), (31) and (33).

a. Best classical stochastic approximation of long-time dynamics generator
We now prove that W defined in Eq. (28) is the closest classical stochastic generator toW in Eq. (26) with respect to the matrix norm induced by the L1 vector norm. For any probability-conserving and positive generator W, we have where in the second and last line we used the positivity of (W) kl ≥ 0, for k = l, the third line follows from the triangle inequality, and the fourth line follows from the probability conservation in both generators The vectorp of barycentric coordinates between metastable phases evolves as  (27)]. As the corresponding state belongs to the MM, at any time t its distance from the simplex of metastable phases is bounded by C cl in Eq.
Let us consider time t such that where c is a constant, so that c √ C cl 1. In this case, the dynamics in Eq. (E4) can be approximated to the linear order with corrections In particular, for the system initially in ρ l that projects ontoρ l in Eq. (14), i.e.,p k (0) = δ kl , from Eq. (E7) we have (E8) On the other hand, from the definition of the absolute value and from the monotonicity of θ(x) ≡ max(x, 0) we have For the classical stochastic generator W defined in Eq. (28), we have from Eq. (E6) that and for the choice c = 1 we obtain Eq. (29).

c. Derivation of Eq. (31)
We have which gives In the second inequality we used the fact that W generates positive trace preserving dynamics and thus e tW 1 = 1, whileW transforms any probability vectors not further away than C cl from the simplex [cf. Eq. (19)], and thus e tW 1 ≤ 1 + C cl . The last inequality follows from Eq. (29). In the first line we used the triangle inequality. In the second line we used that by definition of the projectionP ss on the stationary probability distributionP ss p ss =p ss , and further exploited e tW p ss = p ss and the definition of the induced norm. The last inequality follows from Eq. (31). Note that we did not assume uniqueness of p ss . We note that an related result can be obtained in the non-Hermitian perturbation theory [80], where in the first orderp whereR is the pseudoinverse ofW,RW =WR = I −P ss , with I being the identity matrix (see also Appendix B). Therefore, the first-order corrections can be bounded as [cf. Eq. (29)] for t τ . For the estimate on R 1 2t, where t 2τ , see Eq. (E29) below. Therefore, the condition on the relaxation time in Eq. (32) guarantees that the first-order corrections are negligible.

e. Approximation of dynamics resolvent
Resolvent. We now consider the resolventR of the dynamicsW at 0, i.e., the pseudoinverse ofW,RW = WR = I −P ss , where I is the identity matrix (see also Appendix B). In Eqs. (E22) and (E23) below, we show it can be approximated by the resolvent R of the classical dynamics W in Eq. (28).
We haveR and thus where in the second inequality we used the fact e tWP ss = P ss e tW =P ss so that (e tW −P ss ) ∞ 0 dt (e t W −P ss ) = ∞ t dt (e t W −P ss ). Equation (E19) holds analogously for the dynamics with W. Furthermore, from Eqs. (31) and (33) The first inequality follows from the triangle inequality.
In the last inequality we used the fact P ss − P ss 1 = p ss − p ss 1 and Eq. (E15). Therefore, by applying the triangle inequality we arrive at In the last inequality we used Eqs. (E19) and (E20) together with e tW − P ss 1 ≤ e tW −P ss 1 + e tW − e tW 1 + P ss − P ss 1 2( e tW −P ss 1 + 2 √ C cl t W 1 ) and R 1 ≤ R 1 + R − R 1 . Therefore, The above inequality holds for any time t. When time t can be chosen as in Eq. (32), as in the case of the discussion of Eq. (33), the leading corrections in the righthand side of Eq. (E22) are given by the numerator. In this case, the closeness of resolvents additionally requires that time t can be chosen as where the first inequality must ensure that both e tW −P ss 1 1 and t1 e tW −P ss 1 / R 1 [and thus t is in general larger than in Eq. (32)]. We note that in the second inequality the bound by 1/ √ C cl is relevant for t ≤ R 1 , while at longer times W 1 R 1 / √ C cl is a stricter bound.
Discussion of the results in proximity to a first-order phase transition. We note that both Eqs. (32) and (E23) are fulfilled for all perturbations away from a first-order phase transition that lift the degeneracy of stationary states in the same order. Indeed, in this case C cl and W 1 , as well as 1/ R 1 and 1/τ and 1/τ , are of the same order in the perturbation. Thus,τ W 1 and R 1 W 1 are of zero-order, that is, they remain finite in the limit of the perturbation decreasing to 0, while the upper limits in Eqs. (32) and (E23) diverge in this limit.
Norm of the resolvent. In particular, when the relaxation time τ is defined as the rate of exponential bound on the approach to the stationary state where P ss is the projection on the stationary state, we also have an exponential bound within the MM [see Eq. (D45)] In this case, the norm of the resolvent can be bounded as [cf. Eq. (E18) ] Therefore, if m does not scale with C cl (e.g., is independent of the system size), the condition in Eq. (32) is implied by the condition in Eq. (E23). When the relaxation time τ is defined as we have for the resolvent S of L where in the second inequality we consider t τ (for which e tL − P ss 1). This follows from S = where the second inequality holds for t τ (cf. Appendix D 3).
and otherwise, when (e tW ) kl > 0, where in the last inequality we used e tW 1 ≤ 1 + C cl , as the long-time dynamics transforms the MM onto itself. This gives Eq. (E30).
Finally, T t defined in Eqs. (E34) and (E35) is optimal, as for any discrete trace preserving dynamics T we have T 1 = 1 and thus from triangle inequality Derivation of Eq. (E31). We have where in the second inequality we used Eq. (E30), the fact that T t 1 = 1 from the positivity and the tracepreservation, and e ktW 1 ≤ 1 + C cl as e ktL transforms the MM onto itself.
Derivation of Eq. (E33). We have In the second line we used that the projectionP ss on the stationary probability distributionP ss p ss =p ss , and further exploited T t p ss = p ss and the definition of the induced norm. The last inequality follows from Eq. (E31). We did not assume uniqueness of p ss .

Classical statistics of quantum trajectories
Here, we first consider generalized statistics of jump number, where individual jumps are considered (Appendix E 2 a). Second, we discuss classical approximation of the statistics of a homodyne measurement [73] (Appendix E 2 b) and time-integral of system-observables [74] (Appendix E 2 c) in the presence of classical metastability. Third, we derive corrections to the approximation of generators of statistics of jump number [Eq. (40)], timeintegral of homodyne current [Eq. (E48)] and system observables [Eq. (E63)] (Appendix E 2 d). Finally, we prove the relation of the statistics of jumps and systemobservables in quantum trajectories to statistics of classical dynamics for finite (Appendixes E 2 e and E 2 f) and asymptotic time (Appendix E 2 g).

a. Activity in quantum trajectories
Total activity. In Sec. V B, we discussed the approximation in Eq. (40) of the operator L s in Eq. (36) by the classical operator W s in Eq. (37). This led to the approximation of the maximal eigenvalue θ(s) of L s eigenmode and the corresponding eigenmode ρ ss (s) by the maximal eigenvalue and the corresponding eigenmode W s . Here, we discuss the corrections for those approximations.
We consider LP + (e −s − 1)J as a perturbation of LQ, by means of the degenerate non-Hermitian perturbation theory [80], where Q ≡ I − P is the projection on the fast modes of L [cf. Eq. (4)], with a further (higher order) approximation of LP + (e −s − 1)PJ P = PL s P by W s , which is the superoperator corresponding to W s . The degeneracy of zero-eigenspace P of QL is lifted by W s with the maximal eigenvalue θ (1) (s) [cf. Eq. (42)] corresponding to the eigenmode ρ where (s) j = s j encodes bias values for individual jumps.
In the presence of classical metastability and with internal activity dominating the long-time dynamics, the projection of L s on the low-lying modes [cf. Eq. (38)] where (h s ) j ≡ 1 − e −sj and the individual activities . For the derivation, see Eq. (29) and Appendix E 2 d.
Therefore, θ(s) is approximated by the maximal eigenvalue of W hs in Eq. (E41), while the corresponding density matrix where p ss (h s ) is the maximal eigenmode of W hs . To obtain the relevant corrections, we consider the degenerate non-Hermitian perturbation theory with LP + j (e −sj − 1)J j as a perturbation of LQ, and LP + j (e −sj − 1)PJ j P further approximated by W hs , where J j (ρ) ≡ J j ρJ † j , and W hs denotes the superoperator corresponding to W hs .
Thus, the corrections in Eq. (E44) are given by ρ  ss (s). Here, θ (1) (h s ) denotes the maximal eigenvalue of W hs , which approximates θ(s), with the higher order correction given by Metastability and dynamical phase transitions in activity. For the bias s 1 large enough [but small with respect to the gap to fast modes of the dynamics λ R m − λ R m+1 ], the contribution from W in W hs of Eq. (E42) can be neglected. In that case, m low-lying eigenmodes and eigenvalues of L s are approximated as metastable phases ρ l in Eq. (14).
In particular, for a homogeneous bias, s j = s, and approximation in Eq. (56) we have the correction h s SQJ (ρ l ) − h s (M in −μ in l P) + J SQJ (ρ l ) + (M in −μ in l P) + (L s + h sM in )(ρ l )/h s for the eigenmode R l (s) =ρ l + ..., the correction −h 2 s Tr[P l J SQJ (ρ l )] + (W) ll − h s (J −μ in ) ll to the corresponding eigenvalue θ l (s) = (e −s − 1)μ l + ..., and the cor- Here,M in is the superoperator corresponding toμ in .

b. Homodyne current in quantum trajectories
We now discuss the relation of the homodyne measurement statistics to the classical dynamics between metastable phases.
Statistics of homodyne current. We consider emissions of quanta associated with jump occurrence and denote by dB † j (t) the creation operator for the quanta emitted with the action of jump J j at time t. The statistics of the homodyne current measured with dX(t) ≡ j [e iϕj dB † j (t) + e −iϕj dB j (t)]/2 corresponds to the biased "tilted" master equation [73] where ln{Tr e tLr (ρ) } is the cumulant generating function for the integrated value of the measured homodyne current until time t. The asymptotic statistics is then determined by the eigenvalue θ(r) of L r with the largest real part.
Classical tilted generator for homodyne current. in the presence of classical metastability, for the bias parameter r smaller than the separation to the fast eigenmodes of L, the maximal eigenmode of L r in Eq. (E45) can be approximated as a maximal eigenmode of PL r P, which in the basis of metastable phase corresponds to [cf. Eq. (38)] k, l = 1, ..., m. We haveW r =W − rX + r 2 8 , where (X) kl ≡ j Tr[P k (e −iϕj J jρk + e iϕjρ l J † j )]/2.X can be approximated by the diagonal matrix of the observable averages in the metastable phases wherex l ≡ j Tr[(e −iϕj J j + e iϕj J † j )ρ l ]/2, k, l = 1, ..., m. Therefore, The corrections in Eq. (E48) can be bounded in the leading order by 2 √ C cl W 1 + |r| m j e −iϕj J j max C cl /2 + C + [cf. Eq. (40)]. Here, W r encodes the statistics of time-integral of the system observablex of Eq. (E47) in classical trajectories. From Eq. (E48), θ(r) is approximated by the maximal eigenvalue of W r in Eq. (E48) shifted by r 2 /8. Furthermore, the maximal eigenmode ρ ss (r) of L r is approximated by the correspondingW r eigenmode p ss (r) [cf. Eq. (E64)] as The corrections in Eq. (E49) are given by ρ ss (r)] in the second order of the degenerate non-Hermitian perturbation theory with respect to LP − rX as a perturbation of LQ + r 2 /8 and with a further approximation of LP − rPX P by W r . Here, X (ρ) ≡ j (e −iϕj J j ρ + e iϕj ρJ † j )/2, W r denotes the superoperator corresponding to W r , Q ≡ I − P is the projection on the fast modes of L [cf. Eq. (4)], + stands for a pseudoinverse and S = L + , while θ (1) (r) denotes the maximal eigenvalue of W r . The normalization Tr[ρ ss (r)] = 1 can be achieved by considering the additional correction −Tr[ρ (1) Classical cumulants of homodyne current. For times after the initial relaxation such that τ W 1 t W 1 √ C cl , the rate of the average integrated homodyne current can be approximated as [cf. Eq. (43)] The first term in Eq. (43) where the fluctuations of the average homodyne current in classical trajectories we denoted by the rate of the time-integral of fluctuation rates in metastable phases ( is the average of timeintegral of the homodyne current for lth metastable phase in Eq. (E50), i.e., (p) k = δ kl , k, l = 1, ..., m, and X l ≡ Tr{P l X SQ[ρ(0)]}/Tr[P l ρ(0)] is the average contribution to the integrated homodyne current from before the metastable regime conditioned on the metastable phase that the system evolves into. Therefore, fluctuations of the homodyne current stem from classical transitions between metastable phases with differing average homodyne current and fluctuations inside the metastable phases, corrected by the average from before the metastable regime. When time t in Eq. (E52) can be chosen longer than the final relaxation, the asymptotic rate of homodyne current fluctuations  where l is such thatx l is maximal (for negative r) or minimal (for positive r) among metastable phases. The corrections given in the lowest order of the non-Hermitian perturbation theory are rQSX (ρ l ) − r(X − x l P) + X SQX (ρ l )+(X −x l P) + [L−r(X −X )](ρ l )/r, wherẽ X is the superoperator corresponding tox in Eq. (E47). Furthermore, with the corrections given in the lowest order of the non-Hermitian perturbation theory by −r 2 Tr[P l X SQX (ρ l )]+ (W) ll − r(X −x) ll and by −2rTr

c. Time-integrals of system observables in quantum trajectories
We now discuss the relation of the statistics of timeintegrated system observables to the classical dynamics between metastable phases.
Statistics of time-integrated system observables. The statistics of a system observable M time-integrated over quantum trajectory corresponds to the biased "tilted" master operator [69-72, 74] Here, ln{Tr e tL h (ρ) } is the generating function for time-ordered cumulants of the integral of M until time t. The asymptotic statistics is determined by the eigenvalue θ(h) of L h with the largest real part.
Classical tilted generator for time-integrated system observables. In the presence of classical metastability, for the bias parameter h smaller than the separation of the slow eigenmodes of L to the fast ones, the maximal eigenmode can be approximated as a maximal eigenmode of PL h P, which in the basis of metastable phase corresponds to [cf. Eq. (38)] The corrections in Eq. (E64) are given by ρ The maximal eigenvalues of W h denoted as θ (1) (h) approximates θ(h), with the higher order correction θ (2)   Classical cumulants of time-integrated system observables. For times after the initial relaxation such that τ W 1 t W 1 √ C cl , the rate of the average integrated system observable can be approximated as [cf. Eq. (43)] The first term in Eq. (43) where the fluctuations of the average system observable in classical trajectories is the average of time-integral of system observable in classical trajectories for lth metastable phase in Eq. (E65), i.e., (p) k = δ kl , k, l = 1, ..., m, and M l ≡ Tr{P l MSQ[ρ(0)]}/Tr[P l ρ(0)] is the average contribution to the integrated system observable from before the metastable regime conditioned on the metastable phase that the system evolves into. When time t in Eq. (E52) can be chosen longer than the final relaxation, the asymptotic rate of time-ordered fluctuations is approximated as [cf. Eq. (E48)] In this appendix, we give and prove corrections in approximations ofW s ,W h , andW r . We also discuss how the order of approximation is changed for the metastability in classical stochastic dynamics.
Corrections in classical approximation ofW s . ForW s in Eq. (38), we show below that can be related to the jump activity as where c is a real constant chosen to minimize the operator norm and k, l = 1, ..., m. Equation (E77) illustrates that jump operators considered in the activity can lead to transitions between metastable phases, but not at rates higher than the rates in the effective dynamicsW (up to corrections to the positivity of metastable phases and their projections). Therefore, together with Eq. (29) [and analogously J−(W+μ−μ in ) We also show for individual jumps that [cf. Eq. (E77)] which bounds the corrections in Eq. (E42).
Corrections in classical approximation ofW r . ForW r in Eq. (E48), we show that is approximated by the diagonal matrixx of the observable averages in metastable phases in Eq. (E47) as which together with Eq. (29) gives Corrections in classical approximation ofW h . Similarly, forW h in Eq. (E61) we show below that is approximated by the diagonal matrixm of the observable averages in metastable phases in Eq. (E62) as Therefore, together with Eq. (29), which leads to Eq. (E63).
Derivation of Eq. (E77). We want to compareJ in Eq. (E75) with the effective dynamics in the MM [cf. Eq. (26)] Corrections for classical stochastic dynamics. We now show that metastability in dynamics of probability distributions rather than density matrices (see Appendix B) leads to linear order of corrections in Eqs. (E77), (E79) and (E81).

e. Rates of average and fluctuations in quantum trajectories after initial relaxation
Here, we give corrections to the classical approximations of the average and fluctuations rates of jump number after the metastable regime given in Eqs. (43) and (46). We also discuss the cases of homodyne measurement and time-integral of system observables introduced in Appendixes E 2 b and E 2 c. Derivations can be found at the end of this appendix.
Corrections to average rate of jump number. For time after the initial relaxation, τ t, the rate of average jump number is approximated as in Eq. (43) with corrections [cf. Eqs. (45) and (E174)] where an integer n ≥ 1 is chosen so that t/n belongs to the metastable regime, J = j J † j J j max , μ 1 = max 1≤l≤m |μ l |, and Q ≡ I − P is the projection on the fast modes of L. The first correction arises fro the approximation of the dynamics of fast modes by their complete decay, while the second correction stems from replacing the long-time dynamics in Eq. (26) by the classical stochastic dynamics in Eq. (28) [in the second inequality we used Eq. (29) and (31)]. Since W 1 μ 1 , for the second corrections to be negligible in comparison to μ 1 (an upper bound on K cl (t) /t), we require Furthermore, when |K|/( SQ J ) C n MM , the first correction in Eq. (E104) can be neglected. Alternatively, when where SQ /2 is shorter than any time within metastable regime [see Eq. (E129) below], the contributionK/t from before the metastable regime can be neglected [provided that K cl (t) /t is of the order of μ 1 ], and thus Eq. (E105) is a sufficient condition for the approximation in Eq. (43). In particular, when time t can be chosen longer than the final relaxation t τ [cf. Eq. (32)], we obtain the classical approximation of the asymptotic average rate in Eq. (45). where k, l = 1, ..., m, which can be interpreted as fluctuation rate of metastable phases (see Appendix E 2 f), and (K) kl ≡ δ klKl = δ kl Tr(P l J SQρ) Tr(P l ρ) (E109) encodes the contribution from before the metastable regime. In the second inequality we used Eqs. (31) and (E77). The first and second lines of corrections in the first inequality arise from the approximation of the dynamics of fast modes by the full decay and neglecting any constant contribution to the fluctuations, while the rest of corrections originate from replacing the longtime dynamics in Eq. (26) by the classical stochastic approximation in Eq. (28) (and considering the space of metastable phases with L1 norm, or the space of density matrices with trace norm; cf. Appendix D 3). In particular, we can identify that the third to fifth lines correspond to the corrections to the non-Poissonian classical fluctuations of the activity, [ K 2 cl (t) − K cl (t) 2 − K cl (t) ]/t, the sixth lines is the correction to the time-integral of fluctuation rate inside metastable phases, [ K cl (t) + ∆ cl (t) ]/t, while the seventh and eighth line describe corrections to the contribution from before the metastable regime.
When the corrections in Eq. (E107) are negligible in comparison with the corresponding (positive) contributions to the fluctuation rate, we obtain sufficient conditions on the classical approximation in Eq. (46). As we discuss below, the contributions to the fluctuation rate can be bounded in terms of operator norms leading to generalized conditions, which are sufficient for the classical approximation to hold whenever the bounds saturate (for the opposite case, see e.g., Appendix E 2 f). In particular, if those conditions are fulfilled at time t much longer than the final relaxation, the asymptotic fluctuation rate can be approximated in terms of the classical dynamics [see Eq. (50) and cf. Appendix E 2 g]. Finally, if there exists a leading contribution, e.g., non-Poissonian classical fluctuations characteristic of a first-order dynamical phase transition (see Sec. V B 5), the conditions can be relaxed.
The contribution to the time-integral of fluctuation rate inside metastable phases [ K cl (t) + ∆ cl (t) ]/t is bounded by σ 2 1 , while the contribution to the fluctuation rate from before the metastable regime [the second line of Eq. (E107)] is bounded by 4 μ 1 min( Kp 1 , J SQ ) (see derivation below). Therefore, the first and the second lines of corrections are negligible in comparison with those bounds, when [cf. Eq. (E106)]. Furthermore, the corrections to the time-integral of fluctuation rate inside metastable phases and contribution from before the metastable regime in the last two lines are negligible in comparison with the corresponding bounds when as well as [which is the condition for replacingσ 2 by (σ tot ) 2 , whileμ can be replaced byμ tot in the contribution from before the metastable regime as W 1 / μ 1 2]. Note that Eq. (E112) can be viewed as condition on the maximal anti-bunching of fluctuations inside metastable phases. Finally, the corrections in the third to fifth lines of Eq. (46) are the corrections to [ K 2 They are negligible in comparison with this bound when Eq. (E111) is fulfilled together with It can also be shown that which is a stricter condition than Eq. (E111) for t ≥ 4 R 1 . For those times, no additional condition beyond Eq. (E113) is implied, as the corresponding corrections can be shown to be bounded by Corrections to average rate of time-integrated homodyne current. Similarly to Eq. (E104), for the integrated homodyne current we obtain that the corrections in Eq. (E50) are bounded by where an integer n ≥ 1 is chosen so that t/n belongs to the metastable regime, X ≤ j e −iϕj J j max and x 1 = max 1≤l≤m |x l |.
Corrections to fluctuation rate of time-integrated homodyne current. Similarly to Eq. (E107), for the integrated homodyne current we obtain that the corrections in Eq. (E52) are bounded by and X P min( X , X 1 ). In the first inequality, the first line of corrections in Eq. (E116) arises from neglecting the dynamics of fast modes, while the rest of corrections originate from replacing the long-time dynamics by the classical stochastic dynamics. In the second inequality we used Eqs. (31) and (E81).
Corrections to average rate of time-integrated system observables. Analogously to Eq. (E115), for the time-integral of the system observable the corrections in Eq. (E65) are bounded by where an integer n ≥ 1 is chosen so that t/n belongs to the metastable regime, M ≤ M max and m 1 = max 1≤l≤m |m l |.
Corrections to fluctuation rate of time-integrated system observables. Analogously to Eq. (E116) for the timeintegral of a system observable we obtain that the corrections to the classical approximation of time-ordered correlations in Eq. (E67) are bounded by Note that it follows that Eq. (E127) holds for times t τ , that is, also before the stationary regime. For times after the initial relaxation, τ t, let an integer n ≥ 1 be such that t/n τ belongs to the metastable regime.
where we used Eq. (E128) with n such that t/n belongs to the metastable regime. Furthermore, we recognize and Im(λ k ) not of higher order than Re(λ k ) [generally true for classical dynamics]. In the inequality we used e tL P ≤ 1 + C + together with Eq. (E127) for time SQ . Analogously to Eq. (E133) we have Tr[J Se (t−t1)L QJ e t1L Qρ] can be bounded by J 2 SQ e (t−t1)L Q Qρ or J 2 SQ e t1L Q Qρ (used in the first and the second integrals, respectively). We also assumed that n ≥ 2, so that t/2 τ [cf. Eq. (E128)].
whereK is defined in Eq. (E109) and K (l) cl (t) is K cl (t) for the system in the lth metastable phase during metastable regime, P(ρ) = ρ l , k, l = 1, ..., m.  +2Tr[X S 2 QX e tL Pρ] + 2Tr(X SQX SQρ) where the integer n is such that t ≡ t/n belongs to the metastable regime. Thus, [cf. Eq. (E137)] Rates of average and fluctuations of jump number during metastable regime. For a general initial state ρ, the rate of average number of jumps within the metastable regime, τ t τ , is approximated as [cf. Eq. (E104)] wherep l = Tr(P l ρ) [cf. Eq. (17)], andK ≡ −Tr(J SQρ) is the contribution to the number of jumps from before the metastable regime. We have μ 1 = max 1≤l≤m |μ l |, J = j J † j J j max , and J P ≤ J can be replaced by μ 1 + (C cl + 2C + ) J [100]. We can further replaceμ l byμ tot l introducing additional corrections bounded by √ C cl W 1 . Note that for time we can neglect the contribution from the long-time dynamics (non-Poissonian fluctuations contribution to the classical fluctuation rate) (we then obtain Eq. (53)). Alternatively, this contribution can be neglected when min( J P , m μ 1 ) which condition can be viewed as the lower bound on the anti-bunching of fluctuations inside metastable phases.
Rates of average and fluctuations of jump number in metastable phases. For the initial state ρ chosen as the closest state to the projectionρ l of lth metastable phase [Eq. (14)], we have Therefore, if min(m, μ 1 / J P )(C MM + C + ) 1, and there exist time t C + SQ J / μ 1 [cf. Eq. (E146)], e.g., when μ 1 is of the order of J , the corrections are negligible. We thus conclude thatμ l in Eq. (E145) is the activity of lth metastable phase (see also Appendix E 2 h).
For the initial state chosen as the closest state to the projectionρ l of lth metastable phase, we also have The contribution independent of long-time dynamics, (δσ2 ) l ≡σ 2 l − (μ) l [cf. Eq. (E148)], can thus be viewed as the non-Poissonian contribution to fluctuations of lth metastable phase (see also Appendix E 2 h). We note that σ 2 l can be further replaced by (σ tot l ) 2 introducing additional corrections bounded by √ C cl W 1 [cf. Eq. (29)], which requires Eq. (E114).
The corrections in Eq. (E152) to be negligible in comparison with σ 2 1 , imply a lower and upper bounds on time t, the condition mC MM 1, as well as the following condition on the anti-bunching of fluctuations inside metastable phases, A similar condition in Eq. (E150) guarantees that we can neglect the linear-in-time contribution ( non-Poissonian fluctuations in the classical dynamics). We note that, alternatively, for the approximations in Eqs. (E151) and (E152) we can consider an initial state ρ such that P(ρ) =ρ l , in which case the corrections are bounded as in Eqs. (E104) and (E107), respectively.
For the initial state chosen as the closest state to the projectionρ l of lth metastable phase, we further have [cf. Eq. (E151)] as well as [cf. Eq. (E152)] +4C + X P X SQ .
Therefore, provided that corrections in Eq. (E158) and (E159) are negligible, we can conclude thatx l in Eq. (E155) andσ 2 l in Eq. (E157) are the rates of average and fluctuations, respectively, of an integrated homodyne current in lth metastable phase.
Rates of average and fluctuations of time-integrated system observable during metastable regime. Analogously, for the time-ordered integral of a system observable M we have [cf. Eq. (E118)] Therefore, provided that corrections in Eq. (E164) and (E165) are negligible, we can conclude thatm l in Eq. (E161) andδ 2 l in Eq. (E163) are the rates of average and fluctuations, respectively, of a time-integrated system-observable in lth metastable phase.
Derivation of Eqs. (E144) and (E147). We consider time t within the metastable regime, that is, τ t τ , in order to investigate the cumulants in Eqs. (E130) and (E137) with respect to the metastable state of the system. Therefore, within the MM, we can expand the long-time dynamics in the Taylor series, while the contribution from the outside is bounded as in Eq. (E128).
We have from Eq. (E137) and where we used the Taylor series e tL P = P + tLP + ... assuming t τ [note that Qρ can be replaced by 1 in derivation of Eq. (E130)]. We also used Eq. (E127) so that t L MM J P C MM J P or t μ 1 W 1 ≤ mC MM μ 1 from Appendix (D 3). For time in Eq. (E106) [cf. Eq. (E129)], we can further neglect the contribution from before the metastable regime Similarly, for time in Eq. (E110) we have where we replacedJ =p ss byμp ss in the second line introducing corrections Eqs. (28) and (29) where we used μ tot −μ in Classical approximation of asymptotic rate of integrated homodyne current fluctuations.
Classical approximation of asymptotic rate of integrated system observable fluctuations.

h. Multimodal distribution of quantum trajectories
Here, we show that the distribution of continuous measurement results is generally multimodal during the metastable regime. This is a consequence of distinct continuous measurement statistics in metastable phases of the system. The modes correspond to measurement distributions for metastable phases and the probabilities are determined by the decomposition of the metastable state between metastable phases. See also Sec. V B 4.

Conditional distribution of continuous measurement.
We consider the statistics of jump number conditioned on the final system state at time t in a quantum trajectory evolving into the metastable phaseρ l . This condition can be described as obtaining outcome l the POVM in Eq. (D6) performed on the system at time t, which takes place with the probability approximated by [p(t)] l ≡ (e tWp ) l up to corrections 2 √ is the integral of the total activity in the classical trajectories, conditioned on the system found at time t in lth metastable phase, while the second term encodes the contribution from before the metastable regime [cf. Eq. (46)  , 0), the conditional distribution features the average approximated by the average jump number from a metastable phaseμ l t, while the fluctuation rate is constant. Therefore, with probability approximated byp l , and long enough t τ , the activity k(t) ≡ K(t)/t typically takes values close toμ l (with fluctuations decaying inversely in t), and thus can be interpreted as a single mode in the jump number distribution. When mC MM +C cl 1, there exists at least one such probabilityp l . Therefore, up to corrections 2mC MM + 2C cl + C cl from replacing by 0 allp l that do not fulfill Eq. (E190) (cf. Eq. (C1)), we can understand the probability distribution of jump number as multimodal, with the overall fluctuation rate, Eq. (E147), featuring linear in time contribution due to the different averages of the modes. Note that corrections to multimodal distribution are always present, even for ideal classical MMs (C cl ,C cl = 0), since the long-time dynamics connects the metastable phases (C MM > 0).
In Secs. V B 4 and V B 5, we discuss how measurement of activity of system can be used to identify the metastable phases a stochastic system state corresponds to, provided that metastable phases differ in their internal activity which dominates transition rates of the long-time dynamics.
Finally, we note that analogous arguments hold for homodyne current (cf. Appendixes E 2 b and E 2 f). Moreover, for the constant rates of the conditional average and fluctuation rate, we do no longer require the condition in Eq. (E187).

Cumulants
of conditional finite-time statistics.
Let {P l } l denote a POVM, that is, P l is Hermitian, P l ≥ 0 and l P l = 1 so that p l = Tr(P l ρ) is the probability of obtaining an outcome l in the corresponding generalized measurement on the state ρ. The statistics of jump number conditioned on outcome l at time t is given by [cf. Eq. (E121)] Θ l (s, t) ≡ log Tr(P l e tLs ρ), so that the conditional average of the jump number and its square [cf. Eqs. (E122) and (E122)] Tr(P l e tL ρ) , (E192) Tr(P l e tL ρ) .
We are now ready to approximate the conditional statistics in terms of classical dynamics.
Derivation of Eq. (E178). For a POVM operator P l being a combination of low-lying modes, P † (P l ) = P l , we have that the probability of observing outcome l where (v l ) k = Tr(P lρk ), k, l = 1, ..., m, so that v l max ≤ 1 + C + . In particular, choosing P l as POVM in Eq. (D6), we have that it corresponds to a classical measurement, v l ≥ 0, and further where [p(t)] l ≡ (e tWp ) l . Therefore, the corrections to where we used Eq. (E127) and Eq. (E128) together with (E129) (n ≥ 1 is an integer such that t/n belongs to the metastable regime). Furthermore, [cf. Eq. (E147)]  ≈ Tr(P l e tL J SQJ SQρ) Similarly, Choosing P l as POVM in Eq. (D6), we arrive at [cf. Eq. (E200)] where K 2 cl (t) l is given in Eq. (E181). Similarly, cl (t) l is given in Eq. (E184), as well as where ∆ cl (t) l is defined in Eq. (E182). Eqs. (E202)-(E208) Eqs. (E202)-(E208) describe all corrections to Tr(P l e tL ρ) K 2 (t) l . Therefore, we conclude that the relative corrections to Eq. (E180) are given by the sum of those corrections plus the corrections to the conditional average in Eq. (E178) multiplied by 2Tr(P l e tL ρ) K(t) l , all divided by Tr(P l e tL ρ)[ K 2 (t) l − K(t) 2 l ], plus the correction to the probability in Eq. (E196) divided by Tr(P l e tL ρ). When those corrections are negligible, we arrive the classical approximation in Eq. (E180).
we arrive at the first line of Eq. (E185), and further that W +μ − (J +μ in ) 1 ≤ √ C cl W 1 [cf. Eq. (29)], at the second line of Eq. (E185). We also note that the contribution in Eq. (E185)from before the metastable regime can be neglected for long enough time t τ , as the resulting relative correction is bounded by Tr(P l J SQρ)/(μ l t+ mC MM ).
Finally, from (E202) for time within the metastable regime we have 2(tC MM + SQ )(1 + Qρ ) SQ J P J P l max +2 SQ 2 J 2 Qρ P l max +t 2 C MM min( J P 2 P l max , m J 2 1 v l max ) +2tC MM min( J P J |SQ Qρ P l max , m J 1 Kp 1 v l max ) +tC MM min(2 J P J |SQ P l max , m Σ 2 1 v l max ).
For P l in Eq. (D6), we further have ReplacingW +μ by J +μ in further introduces corrections bounded by √ C cl W 1 (2t J 1 + Kp 1 + 1) [cf. Eq. (29)] and we obtain the classical approximation of the conditional square of jump number. Using then already derived Eq. (E185) we arrive at the approximation of the variance, and thus also the fluctuation rate as given in Eq. (E186).
Derivation of Eqs. (E188) and (E189). When the contribution from before the metastable regime in Eq. (E185) is negligible, the relative corrections from neglecting the long-time dynamics are bounded by W /μ l .
Similarly, in the first line of Eq. (E186) we havẽ l , and thus neglecting the long-time dynamics removes any time-dependence from the fluctuation rate, with the neglected terms bounded in the leading order by Eq. (E127)]. Furthermore, the constant fluctuation rate can be further simplified by neglecting the long-time dynamics in the second term, which leads to corrections bounded by W 1 /p l , and in the third term, with corrections bounded by W 1 Kp 1 /p l . Therefore, we obtain the approximation in Eq. (E189) whenever the relative corrections [2mC MM μ 1 + W 1 (1 + Kp 1 )]/(Σ 2p ) l .

Appendix F: Classical hierarchy of metastabilities
Here, we consider a further separation in the real part of the m low-lying eigenvalues of the master operator, which corresponds to existence of the second metastable regime in the long-time dynamics of the system. We prove that for the classical metastability of m low-lying modes, the second metastability is classical as well.

Hierarchy of metastabilities
Apart from −λ R m −λ R m+1 as considered in Sec. II, here we assume also that −λ R m2 −λ R m2+1 with m 2 < m. Let τ 2 ≥ τ be the timescale of the second relaxation, so that at time t τ 2 the system state can be approximated as [cf. Eq. (3)] (usually τ 2 ≈ −1/λ R m2+1 ). If we can choose time so that the decay of the slower modes is negligible τ 2 t τ 2 (usually τ 2 ≈ −1/λ R m2 ), there exist another metastable rgime during which system states appear stationary [cf. Eq. (4)] with P 2 denoting the projection on the second MM. We denote the corrections to the stationarity in Eq. (F2) as C

Hierarchy of classical metastable manifolds
Assuming the first MM to be classical (see Sec. III), we now show that he structure of the second MM is classical corresponding to m 2 metastable phases. For a related discussion, see Ref. [66].
As discussed in Appendix E 1, the dynamics e tW in the first MM can be approximated at any time t, by the discrete positive and trace-preserving dynamics T t [see Eqs. (E30) and (E31)]. Let t τ 2 in Eq. (E30) be the shortest time t within the second metastable regime. Furthermore, let an integer n be such that nt τ 2 , i.e., nt also belongs to the second metastable regime. We then have whereP 2 denotes the projection P 2 on the second MM in the basis of the metastable phases of the first MM and we introducedC (2) MM ≤ mC (2) MM to denote the corrections to the stationarity in L1 norm (see Appendix D 3). The first inequality follows from the triangle inequality, the second from Eqs. (D41) and (E31). From Eq. (F3) we obtain that the discrete evolution with the transfer matrix T t can be approximated by the same operator for all n min(1/C cl , τ 2 /t). Therefore, for such n the corresponding probability distributions are approximately stationary, that is, metastable (see below the discussion of the figures of merit for the stationarity). From Ref. [30] it is known that metastable probability distributions in classical discrete dynamics can be approximated as mixtures of m 2 approximately disjoint probability distributions (classical metastable phases). Let P t be a projection on the MM of T t and C MM be the corresponding corrections to the stationarity (in the classical space of m configurations corresponding to metastable phases) [cf. Eq. (10)]. From Eq. (F2) we then have We now consider m 2 candidate metastable phases for the second MM. Let p 1 , ... p m2 be the classical metastable phases for T t [cf. Eq. (13)], which under the projection P t on the slow modes of T t we denote p l2 ≡ P t p l2 , l 2 = 1, ..., m 2 [cf. Eq. (14)]. Let v k2 be the dual basis to p l2 , i.e., v T k2 p l2 = δ k2l2 , k 2 , l 2 = 1, ..., m 2 , and the corresponding corrections to the classicality [cf. Eqs. (17) and (19)]. We consider m 2 candidate metastable phases for the second MM as the projections of metastable phases in classical dynamics T t [cf. Eq. (14) and see Appendix C 2] We now estimate the corresponding corrections to classicality. Letp (2) denote the barycentric coordinates of P 2 (ρ) in Eq. (F2) in the basis of Eq. (F6) [cf. Eq. (17)]. We can relatep (2) to the barycentric coordinatesp in the basis of metastable phases p l2 of T t , asp (2) k 2 , l 2 = 1, ..., m 2 . This transformation is close to identity where in the second line we used p l2 1 = 1, while in the last line v k2 max ≤ 1 + C cl /2 [cf. Eq. (F5)] and Eq. (F4). Therefore, where we used p 1 1 + C cl + m 2 C cl [it is 1 + C cl for barycentric coordinates of projected probability distributions over m metastable phases, cf. Eq. (F5), but we projectp with p 1 ≤ 1 + C cl , cf. Appendix (D 3)]. Thus, we arrive at We conclude that when the first MM of an open quantum system is classical, the second metastable manifold is classical as well.
Finally, we note that we can consider a different figure of merit to describe metastability within a given time regime τ t τ , namely, the distance to the best stationary approximation to metastable states By considering ρ as the closest state to P[ρ(0)] and the triangle inequality we obtain a bound in terms of the already introduced figure of merits in Eqs. (10) and (11), We now discuss this figure of merit in the case of hierarchy of metastabilities. Let us denote by C the correction Eq. (F11) in the L1 norm of probability vectors for T n t , where τ 2 nt τ 2 . From Eq. (E31) and Eq. (D41) in Appendix D 3, in analogy to Eq. (F3), we obtain that whereC (2) are the corrections as in Eq. (F11), but in L1-norm and for dynamics withW during the second metastable regime [we have andC (2) ≤ mC (2) , where C (2) are the corrections in Eq. (F11) for τ 2 t τ 2 ]. Therefore, the discrete dynamics T n t exhibits metastability for all n min(1/C cl , τ 2 /t). Furthermore, in analogy to Eq. (F4), the distance of the closest probability distribution to given metastable distribution under dynamics with e ntW and T n t during the (second) metastable regime is bounded as In general the norm W 1 is dominated by the fastest transitions between metastable phases, while the final relaxation time τ describes the slowest timescale of the dynamics. As a consequence, in the presence of the second metastable regime, the classical dynamics W defined in Eq. (28) in general does not approximate the long time-dynamics up to and beyond the relaxation time. Furthermore, the conditions in Eqs. (32) and (E23) are in general no longer fulfilled, and thus the stationary state and the resolvent ofW in Eq. (26) are no longer approximated well by the stationary state and the resolvent of the classical stochastic dynamics W. We now discuss the approximations can be modified to take into account the hierarchy of metastable regimes in the system dynamics.
Approximating the second MM. Instead of considering the approximation of the long-time dynamics by the discrete dynamics T t we can consider a generally weaker approximation by the classical stochastic generator W in Eq. (28). This gives instead of Eq. (F3) [cf. Eqs. (31) and (33)]. Therefore, if the relaxation timẽ τ 2 ≥ τ 2 towards the second metastable regime (with respect to e tW −P 2 1 rather than e tL − P 2 ) fulfills [cf. Eq. (32)]τ the long-time dynamics is well approximated after the relaxation into the second metastable regime. Furthermore, for times the second metastability corresponds then to the metastability in the classical dynamics generated by W, with P 2 approximated by the projection P on the low-lying eigenmodes of W as [Eq. (F4)] where C MM denotes the corrections to the stationarity of e tW , and we considered t τ 2 .
Approximating long-time dynamics for t ≥ τ 2 . After the second metastable regime t ≥ τ 2 , when system states are restricted to the smaller second MM in Eq. (F2), the system dynamics P 2 LP 2 [cf. Eq. (6)] in the basis of m 2 metastable phases denoted byW (2) [cf. Eq. (26)] is approximated by classical stochastic dynamics W (2) [cf. Eqs. (28), (29) and (31)], e tW (2) − e tW (2) where C (2) cl is in the approximation of the second MM by the simplex of m 2 metastable phases and we consider the L1-norm in the basis of m 2 metastable phases. Instead, in the basis of m metastable phases of the first MM, we denote the action of W (2) byW 2 and we have (cf. Appendix D 3) Note thatW 2 in Eq. (F28) is defined on the image of P 2 and in general is not a stochastic generator, as it is probability conserving but only approximately positive.
Approximating the stationary state. When the final relaxation timeτ 2 [with respect to e tW (2) −P (2) ss 1 , wherẽ P (2) ss denotes the projection P ss on the stationary state ρ ss in the basis of m 2 metastable phases; we haveτ 2 ≥ τ (cf. Appendix D 3)] fulfills from Eq. (F27), the stationary state is captured by the stationary distribution p Approximating the dynamics resolvent. We now show that the resolvent of the dynamics can be approximated in two steps: for the faster modes by the resolvent of W and for the slower modes by the resolvent of W 2 , see Eqs. (F38), (F39) and (F40) below (cf. Appendix (E 1)). We where in the first line we introduced the projection P ss,2 on the stationary state of W 2 , and the second and third line refer to operators in the basis of m 2 metastable phases of the second MM. For the resolvent R 2 for the dynamics W 2 , we have [cf. Eq. (F31) and Appendix D 3] The third term in Eq. (F31) can be approximated by the classical dynamics W as [cf. Eqs.
The above inequality holds for any choice of times t 1 ≤ t. When time t can be chosen as in Eq. (F29), and t 1 as time in Eq. (F25), the leading corrections in the right-hand side of Eq. (F38) are given by the numerator as e t1W (I −P 2 ) 1 2C (2) MM [58]. In this case, the small corrections to the approximation of the resolvent additionally require that time t is such that both P (2) ss − e tW (2) 1 1 and t P (2) ss − e tW (2) 1 / R 1 1 [and thus t is in general larger than in Eq. (32)], as well as t 2 C (2) cl W (2) 1 R 1 1, which gives the condition [cf. Eq. (E23)] Finally, time t 1 must be such that t 2 1 √ C cl W 1 / R 1 1, which together with Eq. (F25), gives [cf. Eq. (E23)] (F40) Note that t 1 P −P 2 1 / R 1 1 in Eq. (F38) follows from P −P 2 1 1 in Eq. (F26), since in Eq. (F40) we can choose t 1 not larger than the final relaxation time t 1 ≤τ ≈ R 1 (cf. Appendix E 1).

b. Hierarchy of discrete approximations of classical long-time dynamics
We now consider approximation of the long timedynamics by the classical discrete dynamics T t defined in Eqs. (E34) and (E35). The approximation corrections e tW − T t 1 for time t chosen before the final relaxation are dominated by the fastest transitions between m metastable phases. As a consequence, in the presence of the second metastable regime, T t with time before the relaxation to the second MM, t τ 2 , in general does not approximate the long time-dynamics t ≥ τ 2 [cf. Eqs. (E30) and (E31)]. We now discuss how the approximations can be modified to include the hierarchy of timescales.
Approximating the second MM. Let t τ 2 , possibly shorter than the second metastable regime. If the relaxation timeτ 2 ≥ τ 2 towards the second metastable regime (with respect to e tW −P 2 1 rather than e tL −P 2 ) fulfills [cf. Eq. (E32)] the long-time dynamics is well approximated after the relaxation into the second metastable regime. In particular, the projection of the second MM is approximated by the projection P t on the eigenmodes of T t with absolute value of eigenvalues close to 1 as [cf. Eq. (F4)] P 2 − P t 1 ≤ P 2 − T n t 1 + T n t − P t 1 (F42) C (2) MM + nC cl + C MM 1, where C MM denotes the corrections to the stationarity of T n t .
Approximating long-time dynamics for t ≥ τ 2 . After the second metastable regime t ≥ τ 2 , when system states are restricted to the smaller second MM, the system dynamics e tW (2) in the basis of m 2 metastable phases [cf. Eq. (26)] is approximated by classical discrete dynamics T where C (2) cl is in the approximation of the second MM by the simplex of m 2 metastable phases and we consider the L1-norm in the basis of m 2 metastable phases. Similarly, in the basis of m metastable phases of the first MM, where we denote the action of T cl .
Note thatT t,2 is defined on the image ofP 2 only. Furthermore, although it is trace preserving, it is not positive. We note, however, that since T t can approximate the long times dynamics at times nt, where n 1/C cl [cf. Eq. (E31)], for t chosen for long enough (e.g., t τ 2 ) the dynamics inside the second MM can be captured by T n t .
Approximating the stationary state. When the final relaxation timeτ 2 [with respect to e tW (2) −P (2) ss 1 , wherẽ P (2) ss denotes the projection P ss on the stationary state ρ ss in the basis of m 2 metastable phases; we haveτ 2 ≥ τ (cf. Appendix D 3)] fulfills from Eq. (F43), the stationary state is captured by the stationary distribution p In the basis of m metastable phases of the first MM (or in the trace norm), the distance between the corresponding vectors (or matrices) is in the leading order bounded by p where we used U j 1 ≤ 1 + C cl (as U j transforms the MM onto itself), Π 1 = 1, and Eq. (G15). On the other hand, we have that U n is also a dynamical symmetry [cf. Eq. (61)] and thus from Eq. (G15) there exist a permutation with the matrix Π (n) such that and thus from the triangle inequality Π n − Π (n) 1 ≤ U n − Π n 1 + U n − Π (n) 1 (G21) 7(n + 1)C cl Therefore, for n C cl 1 we can identify so that [cf. Eq. (G20)] Furthermore, for N divisible by n using Eq. (G23) we can retrace its proof with U replaced by U n and n replaced by N/n, in order to arrive at [provided that (N/n)C cl 1]. We conclude that in Eq. (G22) we only require n C cl 1 for all prime factors n of n (and not necessarily n C cl 1).

b. No nontrivial continuous symmetry of metastable phases
As any dynamical symmetry corresponds to an approximate permutation of metastable phases, we argue that continuous dynamical symmetries restricted to low-lying modes are necessarily trivial.
Let us now consider a continuous symmetry, i.e., U φ = e iφG , where with G being a Hermitian operator on the system space. (W ) kl = (W ) π(k)π(l) , where π is the permutation corresponding to Π and k, l = 1, ..., m. We now prove that this condition holds for W as well, as given by Eq. (71). where we introduced k ≡ π(k). This ends the proof [cf. Eq. (B14)].

Example of classicality test with dynamical symmetry
Here, we provide a simple example, why, even in the case of the symmetry fully determining the eigenmodes of the dynamics [cf. Eq. (73)], the classicality test is needed.
We consider a finite system with m disjoint stationary states ρ 1 , ...ρ m and the corresponding projections P 1 ,...,P m [Tr(P k ρ l ) = δ kl , m l=1 P l = 1]. We assume that the stationary states are connected by a dynamical symmetry as U(ρ l ) = ρ l+1 , l = 1, ..., m (with periodic boundary conditions on the label, m + 1 ≡ 1). We have that eigenmatrices of the symmetry fulfill e iϕ k R k = (1/m) m l=1 (e −i2π j k m ) l ρ l and e −iϕ k L k = m l=1 (e i2π j k m ) l P l [cf. Eq. (73) and the normalization L k max = 1], where e iϕ k is an arbitrary global phase, k = 1, ..., m, and e i2π j k m is the corresponding symmetry U eigenvalue, j k ∈ {0, ..., m − 1}. As the invariant set of candidate phases, let us chooseρ (achieved for any of ρ l , l = 1, ..., m), which diverges to ∞ as p → 0 and the basis in Eq. (G56) stops being linearly independent. We note, however, that when expressed in the basis ofρ l , R k is still proportional (with the factor √ m/p for k ≥ 2) to Eq. (73a), as so is L k to Eq. (73b) when expressed in the basis ofP l (with the factor p/ √ m for k ≥ 2).
In this example, the proportionality factor between R k and Eq. (73a) [or L k and Eq. (73b)] indicates that the choice of the basis in Eq. (G56) is not optimal. In the presence of classical metastability, however, rather than for a classical phase transition as considered in this example, the proportionality factor for correctly identified phases does not equal 1, asρ l in Eq. (14) are only approximately disjoint (see Sec. IV), as well asP l are not bounded by 1 or orthogonal (see Appendix D 1). We would expect, however, that the proportionality factor can be related to the corrections in Eq. (19), and should be close to 1 for classical metastable phases.

Rotations of the basis of eigenmodes to expose metastable phases
We now explain how rotations of the basis of eigenmodes allow to solve problem of degeneracy in the coef-ficients, as well as they ensure that a given metastable phases corresponds to extreme value of one of the coefficients when the volume of the simplex of metastable phases in the coefficient space is nonnegligible.
The required lack of degeneracy (up to order C cl ), between metastable phases with the maximal (minimal) kth coefficient in Eq. (H1), corresponds to kth axis in the space of coefficients being normal to the supporting hyperplane at lth vertex of the simplex in the coefficient space where l = l max k (l = l min k ). In convex geometry it is known that for any vertex in a simplex, there exists a supporting hyperplane, that is, there exist a rotation such that the rotated kth axis is perpendicular to a supporting hyperplane. In our case, the lack of degeneracy up to order C cl , additionally requires the distance from that hyperplane of other vertices δ (l) k C cl for l = l max k (l = l min k ). As we argue below, this condition can be translated into a nonnegligible volume of the metastable phases simplex S in the space of coefficients, which is guaranteed by where s m−1 is the maximal volume of a simplex with m− 1 vertices inside the (m − 1)-dimensional unit hypercube.
The condition in Eq. (H3) also guarantees that the separation between any two metastable phases in the space of coefficients equipped with L2 norm is C cl , and thus the metastable phases are distinguishable in the space of coefficients [cf. Eq. (H1)].  (l = 1) orthogonal to the (m−2) dimensional subspace spanned by other vertices,c l with l = 1, l. The distance of other vertices to the supporting hyperplane atc l normal to δc l is equal to L2 norm δc l 2 . On the other hand, the simplex volume Vol(S) = δc l 2 Vol(S l )/(m − 1), where S l is the simplex of all vertices but lth one. Indeed, we have that Vol(S) = det √ C T C/(m − 1)!, where (C) k−1,l−1 = (c l ) k with k, l = 2, ..., m [65], and Vol(S l ) = det (C l ) TC l /(m − 2)!, whereC l is obtained fromC by removing lth column (we assumed l = 1). Therefore, up to rotations, in order for metastable phases to correspond to extreme eigenstates of dynamics eigenmodes, we require a nonnegligible volume of the MM in the coefficient space. Finally, we note that due to the chosen normalization of L k , the MM is enclosed by a (m − 1)-dimensional unit [cf. Eq. (H3)].
We also note that the distance between two vertices, kth and lth, is bounded from below by δc l in the coordinate system shifted to the kth vertex (earlier k = 1). Therefore, it the distance between any two vertices with respect to L2 norm in the space of coefficients is C cl when Eq. (H3) is fulfilled.
Eq. (H7), while the extreme values of the coefficient c I k for L I k can be approximated by considering the sine instead of the cosine. Therefore, degeneracy in the extreme coefficients c R k or c I k can arise in two ways: A. the coefficients c (l) k together with choice of ϕ k lead to degeneracy of extreme values attained by different plane waves, B. some of the cycles contributing to L k are of length longer than n k , so that they have d l /n k -fold degeneracy in the amplitude of the corresponding plane wave [cf. Eq. (73b)].
Case B degeneracy is a direct consequence of from the presence of the dynamical symmetry, as discussed in the first paragraph, and is always present, e.g., in symmetric eigenmodes. Nevertheless, both degeneracy cases can be remedied; see below.
c. Refined metastable phase construction We now explain how to remedy degeneracy of coefficients arising in the presence of a dynamical symmetry in the approach of Sec. VII A. Although any degeneracy can be generally resolved by random rotations as argued in Appendix H 2, here, we explain how the structure of the MM imposed by the dynamical symmetry in the coefficient space can be exploited to simplify the approach; see also the discussion in Sec. VII A.
Case A degeneracy. In general, cycles of various lengths contribute to L k in Eq. (H6). If we consider L k such that there does not exists L l with n l divisible by n k , however, L k is supported on only cycles with the length n k (as there are no longer cycles with the length divisible by n k ). In particular, the number of cycles of length n k equals the degeneracy of the symmetry eigenvalue e inφ k among low-lying eigenmodes. Therefore, the only coefficient degeneracy that can be present belongs to Case A discussed in Appendix H 5 b. We now explain how it can be remedied by an appropriate choice of the phase ϕ k in L R k of Eq. (76). First, let n = 1, 2, ..., n k −1 be such that it is the closest to 1 among e inφ k and with the positive imaginary part. For any ϕ k in Eq. (76), the difference of the maximal and the next in value coefficient for the metastable phases in lth cycle is less than 2 √ 2|c (l) k | sin(δφ k ) 2 , where δφ k = (nφ k mod 2π)/2, but it also can be 0 in the worst case scenario. By considering both ϕ k | be the maximal weight among the plane waves in Eq. (H6). The minimal difference between maximal coefficients for metastable phases in different cycles is at least max[|c (l) k | cos(δφ k /2) − max l =l |c (l ) k |, 0] either for ϕ (1) k or ϕ (2) k . If needed, this bound can be improved by considering additionally ϕ (n) k = δφ k /2 n , n = 1, .., N in Eq. (76), in which case cos(δφ k /2) is replaced by cos(δφ k /2 N +1 ).
Therefore, there is only a single metastable phase corresponding to the maximal coefficient of L R k when the maximal weight |c (l) k | is nondegenerate up to C cl [also when multiplied by cos(δφ k /2)] and when sin(δφ k /2) sin(δφ k ) C cl . The former condition can be achieved for any plane-wave [any l in Eq. (H6)] by a rotation of all L k with the same symmetry eigenvalue provided that the simplex of the coefficients |c (l) k | for those eigenmodes has a nonnegligible volume (cf. Appendix H 2).
Finally, for each cycle, the metastable phase chosen from an extreme eigenstate of (rotated) L R k can be used to recover other elements of the cycle by applying the symmetry U n k times. Therefore, there is no need to consider L I k or other modes supported on those cycles, i.e., L k with a different symmetry eigenvalue e iφ k = e iφ k but with n k = n k . Furthermore, this choice corresponds to the symmetric set of metastable phases in Eq. (66) (cf. Appendix H 5 a).
Case B degeneracy. After finding the extreme states of L k as described in the above paragraph, we still need to consider L l with the symmetry eigenvalue e iφ l and n l such that there exists L k with n k divisible by n l . In that case L l features Case B degeneracy discussed in Appendix H 5 b with respect to metastable phases connected by π n k /n l in the already found cycles [cf. Eq. (H6)].
We are interested in whether the metastable states on which L l is supported via plane waves have been already found by considering L k . This is the case when the degeneracy of e iφ l among the low-lying eigenmodes equals the number of already considered cycles with their length divisible by n l (which is equal the sum of degeneracies of e iφ k for already considered L k such that n k /n l is an integer; without repetitions, i.e., without e iφ k = e iφ k with n k = n k ). Otherwise, L l is supported on cycles that have not been considered yet. In that case, by choosing L l with n l such that the only L k with n k divisible by n l have already been considered, we obtain that all new cycles on which L l is also supported are exactly of the length n l (and their number is the difference between the degeneracy of e iφ l and the sum of degeneracies of relevant e iφ k ). By considering rotations of all eigenmodes with the symmetry eigenvalue e iφ l , a metastable phase in each cycle with the length n l can be found from extreme eigenstates of rotated eigenmodes, as discussed in Case A above. We note that equal mixtures of already found phases in cycles of length n k connected by π n k /n l will also be found here (due to Case B degeneracy; cf. Appendix H 5 a), but those candidate states will be discarded by choosing the candidates yielding the maximal volume simplex. Finally, other elements of the cycles with the length n l can be recovered by applying the symmetry n l times, while other modes supported on those cycles, i.e., L l with e iφ l = e iφ l but with n l = n l , do not need to be considered.
If not all choices of the eigenmodes are exhausted in the way discussed above, i.e., there exist L j with n j < n l , such that n j divides n l , we again repeat the procedure described in the above paragraph, but with respect to all L l and with n l that are divisible by n j . In particular, L j may be discarded if is supported on already considered cycles, i.e., degeneracy of e iφj is equal to the sum of degeneracies of all relevant e iφ l .