Hierarchical classical metastability in an open quantum East model

For the system with $N=3$ spins, a single gap at $m=4$ is present in the spectrum, so that the metastable manifold can be sampled by plotting the coefficients for random pure initial states as in Fig. 3(a). We observe that the metastable manifold is classical, that is, approximated by a simplex, with coefficients of any metastable state approximated by a probabilistic mixture of the coefficients corresponding to the simplex vertices, which describe states with a single or no excitation [cf. Fig. 1(b)]. Since for $m>4$, as relevant for larger system sizes, such a visual verification of metastable manifold classicality is not possible, we instead turn to the recently proposed approach from Ref. Macieszczak et al. (2020), which we sketch now.

For a set of $m$ candidate states $\rho_{1},....,\rho_{m}$, the corresponding metastable states ${\rho_{\text{ss}}}+\sum_{k=2}^{m}c_{k}^{(l)}R_{k}={\tilde{\rho}}_{l}$, $l=1,...,m$, can be considered as the new physical basis replacing the low-lying modes, so that [cf. Eq. (9)]

Here, $\tilde{p}_{l}=\sum_{k=1}^{m}({\mathbf{C}}^{-1})_{lk}c_{k}$ with $({\mathbf{C}})_{kl}=c_{k}^{(l)}$, are the barycentric coordinates with respect to the simplex of ${\tilde{\rho}}_{1}$, …, ${\tilde{\rho}}_{m}$ in the coefficient space. When the distance of barycentric coordinates from probability distributions is negligible, the metastable state can be approximated as a probabilistic mixture of ${\tilde{\rho}}_{1}$, …, ${\tilde{\rho}}_{m}$. If this is true for any metastable state, the metastable manifold is classical and we refer to ${\tilde{\rho}}_{1}$, …, ${\tilde{\rho}}_{m}$ as metastable phases.

For the range of temperatures and field amplitudes corresponding to presence of gaps in the spectrum of the master operators (cf. Fig. 2), using a version of the algorithm from Ref. Macieszczak et al. (2020) (see Appendix A.3 for details), we found sets of $m$ states for which both the average distance and the maximal distance of barycentric coordinates to probability distributions are negligible; see Fig. 4. In particular, Figure 4(a) shows that the metastable manifold with $m=7$ is well approximated by $7$ metastable phases for broad regime of low temperatures and weak coherent field exactly corresponding to the large separation at $m=7$ in the master operator spectrum - with the parameter values above a certain threshold [shown as black line with grey region below; cf. Fig. 2(a)]. Below the threshold, the metastable manifold instead consists of $m=10$ metastable phases [see Fig. 4(b) with the discussed threshold now shown as white dashed line], except for negligibly small values where the separation at $m=10$ in the spectrum also disappears (cf. Fig. 2). These phases remain metastable also above the threshold, but for a smaller range of values of the field and temperature than in the case of $m=7$, which correspond to the hierarchy of metastabilities, i.e., two gaps in the spectrum of master operator at $m=7$ and $m=10$ [cf. Fig. 2(b)]. These results are in agreement with the perturbation theory results derived in Appendix C, which predict emergence of a classical manifold from a quantum metastable manifold at small $\gamma$ and $\Omega$, but at zero temperature and in the absence of the field indicate that softening the constraint leads to quantum decay of excitations within the quantum metastable manifold (cf. Sec. II.3).

We conclude that, at the chosen soft constrain, the metastable manifolds are classical for low temperatures and weak coherent fields (except for their negligibly small values) and there exists an intermediate parameter region with a hierarchy of metastabilities corresponding to two classical metastable manifolds. We will understand the emergence of hierarchy by studying properties of the metastable phases and their long-time dynamics.

We now discuss the properties of the metastable phases whose probabilistic mixtures approximate the classical manifold of metastable states present in the open quantum East model at low temperatures and weak field [cf. Eq. (10)]. We focus on the parameter regime where there exists a gap in the master operators at $m=10$ [below white dashed line in Fig. 4(a)], which captures all values for which a hierarchy of metastabilities is present, but also a region where a gap at $m=7$ is absent [below white dashed line in Fig. 4(b)].

In Fig. 5(a), we show the spin magnetisation along $z$-axis for the metastable phases. For $m=7$ (first two rows, above the white dashed line), the metastable manifold consist of the state with all spins down (no excitation), and six states with a single spin up (a single excitation). For $m=10$, the manifold additionally contains three states with two excitations at maximally separated sites, i.e., followed by two empty sites [see third row in Fig. 5(a)].

As the probability of a spin up or down seems to decrease with the stronger coherent field, we also confirm (see the insets), that the spins in metastable phases are actually aligned with the rotated eigenbasis, $|u\rangle$ and $|e\rangle$, of the stationary state [see Eq. (4) and Appendix B]. Therefore, the metastable phases with no excitations, single excitation and two excitations can be approximately viewed as $|uuuuuu\rangle$, $|euuuuu\rangle$, $|euueuu\rangle$, respectively, and their translations. We obtained such a structure in the first-order perturbation theory with respect to temperature ($|000000\rangle$, $|100000\rangle$, $|100100\rangle$ with translations; see Appendix C.2), and now we confirm it is the case in the presence of the coherent field.

These pure states, however, are not stable since the presence of an excited spin facilitates dynamics on the spin to its right, in turn facilitating dynamics further along the chain: the metastable states thus feature excitations as much separated as possible, so that the relaxation is as slow as possible. Furthermore, the dynamics facilitated by these excitations cause photon emissions from their right neighbour, resulting in a mixed rather than pure metastable state, i.e., $|e\rangle\!\langle e|\otimes|u\rangle\!\langle u|$ replaced by $|e\rangle\!\langle e|\otimes\rho_{\text{ss,1}}$ (cf. Appendix D and Sec. IV.2). This is confirmed by the purity of the metastable phases in Fig. 5(b), where the phases with a single or double excitation feature a purity slightly below $1$, with a lower purity for the state with more excitations. Furthermore, in the first-order corrections due to temperature, purity is lowered proportionally to $\gamma/\kappa$, and Figure 5(b) suggests it is also the case for the coherent field, with the lowest order contribution scaling with $\Omega^{4}/\kappa^{4}$.

Finally, we note that the pure states are exactly orthogonal, and thus the metastable phases are approximately disjoint, as expected from the general theory Macieszczak et al. (2020). Furthermore, the set metastable phases is invariant under the translation symmetry, which is a consequence of the metastable manifold inheriting the symmetry of the dynamics in Eq. (2) with periodic boundary conditions Minganti et al. (2018); Macieszczak et al. (2020).

The effective generator ${\mathbf{\tilde{W}}}$, pictured in Fig. 6(a) encodes all information needed to predict the evolution of the average system state at long times. It conserves the sum of barycentric coordinates, i.e., $\sum_{l=1}^{m}({\mathbf{\tilde{W}}})_{lk}=0$, which is a consequence of the master operator in Eq. (2) being trace-preserving Macieszczak et al. (2020). Although it does not generate positive dynamics (cf. Appendix E.1), its diagonal elements are negative, while its off-diagonal approximately positive, so that it can be approximated by a classical stochastic generator (cf. Appendix E.2). Importantly, Figure 6(b) confirms that the effective dynamics can be approximated by classical stochastic dynamics across the entire metastable region of the parameter space, with the normalised distance of ${\mathbf{\tilde{W}}}$ to the set of classical stochastic generators much smaller than $1$. In fact, this is a consequence of the classicality of the metastable manifold Macieszczak et al. (2020) we discussed in Sec. III.2.

We also note that in Fig. 6(a), the dynamics features the translation symmetry, i.e., $({\mathbf{\tilde{W}}})_{\pi(l)\pi(k)}=({\mathbf{\tilde{W}}})_{lk}$, where $\pi$ is the permutation that the metastable phases undergo under the translation of spins [cf. Fig. 5(a)]. This symmetry is inherited from the translation symmetry of the open quantum East model with periodic boundary conditions Macieszczak et al. (2020). While it reduces the free parameters of the effective dynamics [to $10$ for $N=6$ and $m=10$], it does not guarantee the presence of detailed balance we demonstrate next.

In the effective dynamics, the stationary probability current between the $k$th and $l$th metastable phase is given by $({\mathbf{\tilde{W}}})_{kl}({\mathbf{\tilde{p}}_{\text{ss}}})_{l}-({\mathbf{\tilde{W}}})_{lk}({\mathbf{\tilde{p}}_{\text{ss}}})_{k}$, where ${\mathbf{\tilde{p}}_{\text{ss}}}$ is the stationary distribution of ${\mathbf{\tilde{W}}}$ (or equivalently the barycentric coordinates for ${\rho_{\text{ss}}}$). Detailed balance is then defined to be when a systems stationary state exhibits no currents (see Appendix E.3).

As a first check of detailed balance in the effective dynamics, we consider the similarity transformation which renders classical detailed-balance generators symmetric,

For the effective generator in Fig. 6(a), we indeed obtain an approximately symmetric matrix in Fig. 7(a).

To verify detailed balance across the range of parameters for which metastability occurs, we consider in Fig. 9(b) the ratio of the total stationary current

to the total activity

which ratio bounds the normalised distance to the closest detailed balance dynamics (see Appendix E.3). We observe that the current across all metastable parameters is small compared to the system’s activity, and thus the long-time dynamics can be well approximated by dynamics with detailed balance. This is also the case for the classical East model ($\Omega=0$), which features only approximate detailed balance when restricted to the metastable manifold [see the bottom panel in Fig. 7(b)], although there are no stationary currents between the $2^{N}$ configurations of up and down spins in the classical system.

For the classical model, these results can be traced back to the perturbative dynamics between configurations. The perturbation effect of the soft constraint removes one excitation at a time with rates proportional to $(1-p)^{2}\kappa$, or reintroduces, removes or shifts a single excitation, at rates proportional to $(1-p)^{2}\gamma$. For the small temperature, phases with double excitations are reduced to a single excitation at rates proportional to $\gamma^{2}/\kappa$, while at rates proportional to $\gamma^{3}/\kappa^{2}$ a second excitation can be introduced or removed, or a single excitation can be shifted. The result is a ladder structure of the dynamics with respect to the number of excitations which necessarily implies detailed balance, though the approximation worsens for larger $\gamma$ or $(1-p)$ due to higher order corrections; see Appendices C.2 and C.3 for details.

Approximate detailed balance observed also in the presence of a weak coherent field in Fig. 7(b), suggests that a similar mechanism may be responsible for the long-time dynamics in the open quantum East model. This is indeed confirmed for the parameters chosen in Fig. 6(a): the most probable transitions (yellow-light green) are associated with the removal of the second excitation, or removal of a single excitation towards the unexcited state; the less likely transitions (green) correspond to a shift of a single excitation, the introduction of one excitation, or removal of two excitations; while the least likely transitions (blue) correspond to the introduction of two excitations simultaneously, or shift of two excitations [see also Figs. 9(b) and 9(c)].

Since the long-time dynamics possesses approximate detailed balance, the stationary state of the effective master operator is effectively thermal. We discuss now its effective free energy function.

In the full classical East model, whose stationary state is a product state, the free energy is simply a function of the number of excitations, but not their relative distance (i.e., it is not dependent on any type of interactions). Thanks to the exponential form of a thermal distribution, this can be tested by considering ratios of state probabilities: the exponents only state dependence will be a linear function of the number of excitations. We test this property in Fig. 8(a) for the distribution over the metastable phases of the stationary state of the the open quantum East model, by comparing the ratio of probabilities of finding the stationary state in one of the single or double excited states, $\tilde{p}_{2}/\tilde{p}_{1}$, to the ratio of finding it in the unexcited state or one of the single excitation states, $\tilde{p}_{1}/\tilde{p}_{0}$. We would expect these ratios to differ when interactions contribute to the effective energy of the phase, due to the presence of multiple excitations in the phase with probability $\tilde{p}_{2}$. However, we find that the ratio of these ratios is close to $1$ at all metastable parameters, suggesting that in the quantum model interactions do not play a role in the free energy of the metastable phases, as in the classical East model.

Furthermore, the free energy per excitation is directly determined by the ratio of probabilities for excited and unexcited states in the single-site dynamics; see Fig. 8(b). This is again the consequence of the product structure of the stationary state, and the metastable states being approximated simply as pure states with a single or double excitations [cf. Fig. 5(a)], which appear as the leading order corrections to the stationary state above the no-excitation contribution ¹¹1In the case of the classical model ($\Omega=0$), in the perturbative regime, the metastable manifold is known to consist of isolated excitations at any system size (see Appendix C), and thus the free energy over metastable phases is again only the function of number of excitations, and so we can expect, also in the open quantum East model for larger system sizes. Therefore, the error of this non-interacting approximation of the free energy will increase as temperatures and coherent field values grow [cf. Fig. 8(b)].

The effective lifetimes of the metastable phases, i.e., the distinct periods in quantum trajectories, are given by the inverse magnitude of the diagonal elements of the effective generator. From the translation symmetry of the metastable manifold, the lifetimes of metastable phases connected under the symmetry must be the same, i.e., we have: the lifetime $\tau_{0}$ of the homogeneous unexcited phase, the lifetime $\tau_{1}$ of six phases with a single excitation, and the lifetime $\tau_{2}$ of three phases with double excitations [see Fig. 9(a)].

The unexcited state is the longest lived metastable phase, as it can only be excited via to the softness of constraint in $H_{j}$ or $J_{j}^{+}$ (for any $j=1,...,N$) [cf. Eq. (6)], which takes place at the rate proportional to $N(1-p)^{2}(\Omega^{2}/\kappa+\gamma)+...$ in the perturbative regime $|1-p|\ll 1$ (see Appendix C.3.2). The simple dependence of the rate on the system size $N$ is due to the translation symmetry: the excitation can be introduced at any of the $N$ spins. In contrast, the removal of a single excitation of $j$th spin, again requires the softness of constraint in $J_{j}^{-}$ (or in $H_{j}$, which however is of lower order), and thus takes place at the faster rate $\kappa(1-p)^{2}+...$, leading to the separation of timescales $\tau_{1}/\tau_{0}\approx N(\Omega^{2}/\kappa^{2}+\gamma/\kappa)$. Furthermore, when two gaps are present in the master operator spectrum [above the white dashed line in Fig. 9(a)], the hierarchy of metastabilities ($m=7,10$; cf. Sec. III.2) is manifested in the distinct values of the lifetimes $\tau_{1}$ and $\tau_{2}$, while for a single gap ($m=10$; below the white dashed line), these lifetimes are necessarily comparable, which we now explain.

The softness of constraint causes the removal rate of an excitation from metastable phases with a single excitation and double excitation to be the same (except from the fact that two possible sites to decay in the double excited state), while for hard constraint only the second excitation can be removed due to the temperature of the coherent field (by flipping the unexcited spins between two excitations which allows for the hard constraint to be fulfilled). When the constraint is soft, the absence or presence of separation between $\tau_{1}$ and $\tau_{2}$ is determined by the softness-induced and temperature-induced dynamics being faster, respectively. In Fig. 9(a), the regime of smaller (greater) $\gamma$ and $\Omega$ below (above) the threshold corresponds to the former (latter) process being the fundamental mechanism in relaxation of the metastable phase with double excitation towards the stationary state.

In Figs. 9(b) and 9(c), we show examples of the effective master operator for two sets of parameters, indicated by the crosses in Fig. 9(a), corresponding to the cases with a single metastability and a hierarchy of two metastabilities. In both cases, a double-excitation phase is most likely transformed into one of two single-excitation phases (equally likely due to the translation symmetry by $3$ sites of the double-excitation phase), with one excitation inducing relaxation of the other. A single excitation phase is most likely transformed into no-excitation phase, which in turn gets excited most likely with only a single excitation (into one of six single-excitation phases). This ladder structure of the effective classical dynamics supports detailed balance in the dynamics discussed in Sec. IV.1.

Beyond those leading order transformations, however, shifts of a single excitations are possible due to two different mechanisms. In Fig. 9(b), we observe that the shift of a single excitation is possible to all sites except that corresponding to the possible position of a second excitation, in which case the second excitation is introduced instead. This indicates that the shift is actually facilitated by the introduction of excitations and their subsequent decay, which can be facilitated either by several excitations by the temperature/coherent field, or the softness of constraint allowing for introduction of excitations directly in unexcited sites. The uniform probability of shifts to different sites in Fig. 9(c), suggest that the two processes contribute equally, while, in the case of hierarchy in Fig. 9(b), the larger values of the coherent field and temperature dominate the latter process and only some shifts are possible. This is directly supported by the perturbation results in the classical model; see Appendix C. We note however, that for considered system size of $N=6$ and the chosen constraint with $p=0.99$, we do not yet capture the hallmark behaviour of the classical East model, where required order of temperature contributing to the dynamics of excitations scales logarithmically with their distance (cf. Appendix C.2.2). In particular, we cannot verify whether the necessarily (quadratically) higher orders in which the local coherent field contributes to the dynamics of in the open quantum East model (cf. Appendix C.4.2) alter this characteristic.

Although there is no apparent directionality in the dynamics for both cases, which is likely due to the small system size $N=6$ (cf. Appendix C.2.4), we would in general expect this to follow from the presence of the constraint to the left [cf. Eq. (6)], and such directionality is present in perturbation theory with respect to temperature for larger system sizes. Nevertheless, we observe in both cases that the unexcited metastable lifetime is much longer than the metastable phases with a single or double excitations; in trajectories of the classical effective dynamics most time is spent in the unexcited state, and excitations are present at isolated moments in time. These periods are also isolated in space, due to the symmetry structure of the metastable manifold (see also Appendix D).

We now discuss how metastability and the structure of long-time dynamics manifests itself in the emission patterns in individual experimental realisations of the system dynamics Olmos et al. (2012); Lesanovsky et al. (2013).

Consider first dynamics in the case of the state being (on average) in one of the metastable phases featuring a single or double excitations. An excitation of site $(j-1)$ fulfils the hard constraint of the single spin dynamics of the $j$th spin [cf. Eqs. (3) and (6)], enabling dynamics on this site and thus its relaxation towards the single-spin stationary state in Eq. (4). Thus, for times shorter than relaxation of the considered metastable phase, $t\lesssim\tau_{1},\tau_{2}$, in an individual realisation of an experiment (or quantum trajectories obtained in QJMC simulations) photon emissions occur from $j$th spin corresponding to the jump $J_{j}^{-}$ [Eq. (6b)], so that the metastable phase with an excitations appears locally bright. These emissions occur at the rate $\kappa{\mathrm{Tr}}(S^{+}S^{-}\rho_{\text{ss,1}})\approx\Omega^{2}/\kappa+\gamma+...$, so that the site next to the excitation appears brighter for higher temperature or coherent field values [see Figs. 9(d) and 9(e)]. In contrast, for the unexcited metastable phase, the hard constraint in the dynamics is not fulfilled, and therefore, this phase appears dark in quantum trajectories, before it relaxes due to the soft constraint introducing of a single excitation at $t\gtrsim\tau_{0}$.

At longer times, higher order processes introducing several excitations or exploiting softness of constraint become non-negligible on average, contributing to the long-time dynamics of the metastable phases by connecting disjoint parts of state space. In individual quantum trajectories these processes take place separately and at fluctuating times, but are typically followed by the fast decay of excitations due to satisfied hard constraints (on timescales $t\lesssim-1/\lambda_{m+1}^{R}$; $m=10$) towards another metastable phase. Therefore, a time coarse-graining of quantum trajectories leads to the system transitioning only between metastable phases. As averaging over trajectories returns the evolution with the master operator [Eq. (1)], transitions in coarse-grained quantum trajectories must be governed by the effective long-time generator [Eq. (13)]. This is corroborated in Fig. 9(d) where, correspondingly with the transition rates of the effective stochastic generator, in the case of a single metastable regime, mostly transitions between the excited states and dark states are observed. Meanwhile, in Fig. 9(e), for the hierarchy of two metastabilities, there is also a significant presence of transitions between states with a single excitation, shifting the location of emissions.

We conclude that the dynamical heterogeneity in the quantum trajectories is the microscopic counterpart to the classical stochastic jumps between phases with different numbers of excitations at different sites, which arise as a result of temporal coarse-graining of quantum trajectories. This is a general phenomenon, detailed theoretically in Ref. Macieszczak et al. (2020). In particular, this relation could be used to explore possible differences in the contributions to the dynamics from the temperature and the coherent field at larger system sizes accessible in QJMC simulations.

Systems with intermittent dynamics are commonly found to exist near a dynamical phase transition in the statistics of the activity, i.e., the number of jumps per unit time Garrahan and Lesanovsky (2010); Garrahan et al. (2011); Ates et al. (2012); Lesanovsky et al. (2013). Here, we demonstrate this for the global activity for jumps related to loss of excitations. Since our system exhibits dynamical heterogeneity, we also find the system in proximity to transitions in the statistics of the local activity.

Dynamical phase transitions in global jump activity. The intermittent emissions in trajectories have a direct effect on the time integrated statistics of their corresponding jumps. The statistics of a trajectory-observable chosen as the number $K^{-}(t)$ of $J_{j}^{-}$ jumps up to time $t$ summed across all sites, is encoded by the cumulant generating function

where

can be obtained using the biased master operator

[cf. Eq. (2)]. Furthermore, the statistics of the activity $k^{-}(t)=K^{-}(t)/t$ is encoded in the long time limit by the scaled cumulant generating function (SCGF)

given by the leading eigenvalue of $\mathcal{L}_{s}$ [cf. Eqs. (17) and (18)]. The corresponding eigenmode ${\rho_{\text{ss}}}(s)$ of $\mathcal{L}_{s}$ is the average asymptotic state in the biased trajectory ensemble, where each trajectory is weighted by $e^{-sK^{-}(t)}$, before the overall ensemble is then renormalised [cf. Eq. (18)]. The SCGF plays the role of free energy in non-equilibrium statistical mechanics Touchette (2009), with its non-analyticities corresponding to dynamical phase transitions Garrahan and Lesanovsky (2010); Garrahan et al. (2011); Ates et al. (2012); Lesanovsky et al. (2013).

In Fig. 10(a), two sharp changes are found in the first derivative of $\theta(s)$, i.e., the average activity $k(s)=-\mathrm{d}\theta(s)/\mathrm{d}s$, at negative bias $s$ close to $0$, between the values equal zero, one or twice the average single spin activity proportional to $\Omega^{2}/\kappa+2\gamma+...$ [see Eq. (4)]. This indicates that the proximity of the physical dynamics $s=0$ to two first-order dynamical phase transitions. Furthermore, these changes occur as the perturbation due to the bias becomes larger than $\lambda_{m}$ for $m=7$ and $m=10$, which indicates their relation to the presence of the hierarchy of metastabilities. This is also supported by Fig. 10(b), where in a decomposition of ${\rho_{\text{ss}}}(s)$ between metastable phases (in its barycentric coordinates), at $s=0$ it corresponds mostly to the dark metastable phase, while for a large enough negative bias $s$ (towards more active trajectories) can be characterised as the equal mixture of six single-excitation metastable phases ($1/6$th probability each) or the equal mixture of three double-excitation metastable phases ($1/3$rd probability each). This homogeneity follows from the translation symmetry of ${\mathcal{L}}_{s}$.

These results actually follow from the correspondence of $\theta(s)$ to an SCGF for integrated metastable phase activity in classical trajectories of the long-times dynamics Macieszczak et al. (2020) (cf. Ref. Macieszczak et al. (2016b)), which holds for small to moderate values of $s$ [when contributions from fast modes are negligible; cf. the inset in Fig. 10(b)] and metastable phases distinguished by the average jump activities dominating rates of the long-time dynamics. This is exactly the situation in the open quantum East model due to the constraint in $J_{j}^{-}$ fulfilled by excitations present in metastable phases [cf. Fig. 5(a)]. Indeed, biasing trajectories with negative $s$ effectuates post-selection towards more active trajectories, which in this case correspond to metastable phases with a single or double activity for smaller or larger $|s|$ (respectively, in order to make up for the shorter lifetime $\tau_{2}\ll\tau_{1}$ due to the hierarchy of metastabilities present for the chosen parameters; probability of trajectories with even higher activity remains negligible). In contrast, for positive $s$ inactive trajectories are preferred, corresponding to the dark metastable phase with no excitations.

Dynamical phase transitions in local jump activity. Rather than the global activity of jumps across the system, we can consider local jump activity, with the locally biased master operator [cf. Eq. (19)]

encoding the joint statistics of number $K_{j}^{-}$ of jumps $J_{j}^{-}$ at sites $j=1,...,N$ observed up to time $t$. Similarly to Fig. 10 for the full jump statistics, in Fig. 11, we observe sharp changes in the first derivative of the corresponding maximal eigenvalue of ${\mathcal{L}}_{s_{1},...,s_{N}}$.

In Figs. 11(a) and 11(b), we look at a cross-section with $s_{1}=s$ and $s_{j}=0$ for $j\neq 1$. This biases towards trajectories containing significant periods of the single excitation metastable phases, and ignores the double excitation phases: indeed, in Fig. 11(a) there is only a single jump in the activity to a value corresponding to the activity of single excitation phases, while the overlap with the phase featuring an excitation at site $6$ that induces emissions on site $1$, turns out to be dominant at negative values of $s$ in Fig. 11(b).

To target a double excitation state, we look in Figs. 11(c) and 11(d) at a cross-section with $s_{1}=s_{4}=s$ and $s_{j}=0$ for $j\neq 1,4$, targeting the phase with excitations on sites $3$ and $6$. As expected, there is a pair of jumps in the activity in Fig. 11(c), corresponding to the activity of single excitation phases and double excitation phases respectively. For smaller negative values of $s$, the overlap with the metastable phases in Fig. 11(b) is split evenly across the single excitation phases on sites $3$ and $6$, as expected in comparison with Fig. 10; for large values the only relevant overlap becomes the double excitation phase that was targeted with this choice of bias.

Metastable phases from biased trajectories. Beyond clarifying the relation of first-order dynamical phase transitions to metastability, Figs. 10 and 11 demonstrate that metastable phases differing in activity can be obtained as the asymptotic average states in the biased ensemble of trajectories. This result indicates an alternative method to uncover the structure of metastable manifold, which does not require the diagonalisation of the master operator (cf. Appendix A.3). While methods for efficient sampling of biased classical trajectory ensembles are common Crooks and Chandler (2001); Bolhuis et al. (2002); Lecomte and Tailleur (2007); Giardina et al. (2011); Ferré and Touchette (2018); Rose et al. (2020), with some work in this direction for quantum systems Schile and Limmer (2018); Carollo and Pérez-Espigares (2020), more development is needed to achieve the speed needed for many-body quantum systems. A possible direction could be the use of tensor network techniques, as done in recent classical large deviation studies Bañuls and Garrahan (2019); Causer et al. (2020).