Quantum jump Monte Carlo approach simpliﬁed: Abelian symmetries

We consider Markovian dynamics of a ﬁnitely dimensional open quantum system featuring a weak unitary symmetry, i.e., when the action of a unitary symmetry on the space of density matrices commutes with the master operator governing the dynamics. We show how to encode the weak symmetry in quantum stochastic dynamics of the system by constructing a weakly symmetric representation of the master operator: a symmetric Hamiltonian, and jump operators connecting only the symmetry eigenspaces with a ﬁxed eigenvalue ratio. In turn, this representation simpliﬁes both the construction of the master operator as well as quantum jump Monte Carlo simulations, where, for a symmetric initial state, stochastic trajectories of the system state are supported within a single symmetry eigenspace at a time, which is changed only by the action of an asymmetric jump operator. Our results generalize directly to the case of multiple Abelian weak symmetries.


I. INTRODUCTION
Markovian open quantum systems describe a broad class of systems interacting weakly with environments whose dynamics are much faster than those of the system itself, as relevant, e.g., for atomic, molecular, and optical physics [1], as well as optomechanics [2]. This leads to system dynamics efficiently described by a local-in-time master equations [3,4], so that both the dynamics and stationary states can be found by its numerical integration or diagonalization. Since the space on which the corresponding master operator acts scales quadratically with the system dimension, other methods for exact numerical simulations of dynamics have been developed scaling with respect to the system dimension rather than its square, such as the quantum jump Monte Carlo (QJMC) approach [5][6][7][8][9], also known as the quantum trajectory technique or Monte Carlo wave-function method, which corresponds to the stochastic description of system dynamics in terms of quantum Langevin equations [10][11][12] or continuous measurement theory [13,14].
Similarly as in closed quantum system dynamics governed by Hamiltonians, the presence of symmetries in master equations is known to simplify the structure of corresponding master operators, although, due to the presence of dissipation, their symmetries are not in general related to conservation laws [15][16][17][18] and their stationary states are unique [19][20][21][22]. In this work we show how a weak symmetry of the master equation can be encoded in the corresponding stochastic dynamics of an open quantum system: via a symmetric Hamiltonian and jump operators connecting only the symmetry eigenspaces with a fixed eigenvalue ratio, which we refer as a weakly symmetric representation of a master operator. This has direct consequences for the numerics: QJMC simulations are simplified, particularly for symmetric initial states, which remain symmetric and thus confined to a single symmetry eigenspace at a time. In turn, also the construction of the master operator, which describes change in time of the average system state, is simplified. Our results carry on directly to the case of multiple weak symmetries, provided their action on density matrices commutes, that is, they correspond to an Abelian group. This is illustrated with a dissipative spin system featuring both a weak translation symmetry and a weak rotation symmetry.
This article is structured as follows. In Sec. II we review Markovian dynamics with weak symmetries. We define weakly symmetric representations and show how to construct them in Sec. III. We then explain how such representations simplify stochastic dynamics in Sec. IV, leading to simplified construction of the master operator and QJMC simulations, as outlined in Sec. V. Finally, we discuss examples of manybody systems with weak symmetries in Sec. VI.

II. WEAK UNITARY SYMMETRIES OF OPEN QUANTUM SYSTEM DYNAMICS
Here we review Markovian dynamics of open quantum systems featuring weak symmetries.

A. Open quantum system dynamics
The Markovian dynamics of an open quantum system is governed by a Gorini-Kossakowski-Lindblad-Sudarshan master equation [3,4], and J j are so-called jump operators describing the interaction with external environments. The dynamics in Eq. (1) is completely positive and trace preserving [L † (1) = 0], from which it follows that there exists a stationary state of the system [L(ρ ss ) = 0]. We refer to the superoperator L as the master operator. A Hamiltonian and jump operators are not uniquely defined for a given master operator [23], and we refer to their particular choice as its representation.

Definition
Dynamics of an open quantum system features a dynamical symmetry [15] or a weak symmetry of the dynamics [16,17] when it commutes with a unitary transformation of system states. For a unitary operator U on the system Hilbert space H, the dynamics features the corresponding weak symmetry when the master operator is symmetric, with respect to the action of the symmetry on density matrices, This is often referred to as a discrete weak symmetry, as it follows that the dynamics features weak symmetries for unitary operators U n , with n ∈ Z. Abelian weak symmetries correspond to weak symmetries with commuting symmetry superoperators, [U 1 , U 2 ] = 0, which requires the symmetry operators to commute as well, [U 1 , U 2 ] = 0. A special case is a continuous weak symmetry, that is, a weak symmetry for U φ = e iφS , where S is a Hermitian operator and φ ∈ R, which requires [cf. Eq. (2)] Note that weak symmetries in Eqs. (2) and (3) are defined as symmetries of the master operator in Eq. (1). We will discuss the resulting properties of a Hamiltonian and jump operators in Sec. III, and the structure of the corresponding quantum stochastic dynamics in Sec. IV.

Resulting structure of master operator
It is known that weak symmetries limit the structure of the master operator, which can be exploited to simplify its diagonalization or numerical integration required to solve the dynamics of system states [17], as we review below.
The structure of the master operator due to a weak symmetry can be conveniently seen in the Liouville representation (see Fig. 1 FIG. 1. Dynamics with a weak symmetry. (a) Consider a basis of eigenstates of a unitary operator U (or a Hermitian operator S) and denote H k the eigenspace with e iφ k (or s k ) eigenvalue. A density matrix ρ features both support inside H k (solid diagonal blocks) and coherences between symmetry eigenspaces (shaded off-diagonal blocks). (b) When the density matrix is expressed as a vector ρ, the master operator becomes a matrix L [cf. Eq. (5)]. The weak symmetry implies that the eigenspaces of the symmetry superoperator U (or S) evolve independently [gray blocks (with black solid or black dashed contour)]. When ratios of U eigenvalues (or gaps in S eigenvalues) are trivially degenerate [e i(φ k −φ l ) = e i(φ k −φ l ) (or s k − s l = s k − s l ) for k = l implies k = k and l = l ], individual coherences evolve independently (gray blocks with black dashed contour vanish while colored blocks remain).
is governed by the matrix where * and T denote the complex conjugation and the matrix transposition in the chosen basis, respectively, and we used where U ≡ U ⊗ U * . Therefore the weak symmetry is equivalent to conservation of the eigenspaces of symmetry superoperator U, where H k is the eigenspaces of U corresponding to an eigenvalue e iφ k , and 1 H k is the corresponding orthogonal projection, while e iδ is an eigenvalue of U. In particular, considering e iδ = 1 leads to k = l in Eqs. (7), we obtain that the symmetric part of a system state evolves independently from the coherences between symmetry eigenspaces.
As a result, in a generic case when the stationary state of the dynamics is unique [19][20][21], it is symmetric and can be found by solving the dynamics restricted to the symmetric eigenspace [17,22,24,25]. Furthermore, averages and higherorder correlations for symmetric system observables can be found, without loss of generality, by solving the dynamics restricted to that eigenspace, that is, by considering symmetric initial states. Finally, choosing the basis as an eigenbasis of U , the matrix U becomes diagonal. Thus, after reordering the basis {|ψ k ⊗ |ψ l } dim(H) k,l=1 of H ⊗ H to group together elements corresponding to the same eigenvalues of U, L in Eq. (5) becomes block diagonal, with the blocks corresponding to the eigenspaces of U (see Fig. 1). The dynamics of the system states can then be found by diagonalization or numerical integration of individual blocks. 1 Similarly, as superoperators of Abelian symmetries commute in the Liouville representation, [U 1 , U 2 ] = 0, the master operator L featuring corresponding weak symmetries conserves all intersections of their eigenspaces. Let orthogonal subspaces {H k } k be defined as intersections of eigenspaces of the symmetry operators, so that each H k corresponds to a different set of eigenvalues for the symmetry operators. Then the sums in Eq. (7) are replaced by sums over k, l with the ratios of the symmetry operator eigenvalues corresponding to a given set of eigenvalues for the corresponding symmetry superoperators. In particular, for a continuous weak symmetry in Eq. (3) and an eigenspace H k of S corresponding to an eigenvalue s k , e i(φ k −φ l ) = e iδ in Eq. (7) is replaced by s k − s l = δ, where δ is an eigenvalue of S = S ⊗ 1 − 1 ⊗ S * , so that δ = 0 corresponds to the symmetric part of a system state.
As we show in Sec. V A, the Liouville representation of the master operator in the symmetry eigenbasis can be efficiently constructed using a weakly symmetric representation. Even when restricted to symmetric states, however, the master operator acts on the space of dimension k dim(H k ) 2 , which can inhibit its diagonalization or numerical integration. Therefore, in Sec. V B, we instead focus on exploiting Abelian weak symmetries to simplify QJMC simulations.

III. WEAKLY SYMMETRIC REPRESENTATIONS
Here, we define and construct weakly symmetric representations for any master operator with a weak symmetry by modifying its Hamiltonian and jump operators into a form respecting the symmetry. We also discuss their nonuniqueness.

A. Definition
In Sec. III B we show that in the presence of a weak symmetry in Eq. (2), there exists a weakly symmetric representation with a HamiltonianH and jump operators {J j } j such 1 The preservation of operator hermiticity by the dynamics and symmetry superoperators, [L(ρ )] † = L(ρ † ) and [U (ρ )] † = U (ρ † ), leads to the dynamics of the blocks corresponding to complex-conjugate pairs of symmetry eigenspaces being related by the complex conjugation and swap operation T , as T (Lρ ) Since the Hamiltonian and the jump operators in Eq. (8) are eigenmatrices of the symmetry superoperator U , they are supported on the corresponding eigenspaces, for all j. This property plays a crucial role in simplifying numerical simulations of the system dynamics (see Sec. V). Note that when e iδ j = 1 for all j, a Hamiltonian and all jump operators themselves are symmetric. This then holds for any representation of the master operator [23] and is known as a strong symmetry [16,17]. Although a strong symmetry implies the corresponding weak symmetry, the converse is not true [15][16][17], as evident by considering weakly symmetric representations in Eqs. (9). Similarly, for Abelian weak symmetries, a Hamiltonian can be chosen symmetric with respect to all symmetry superoperators and jump operators can be chosen as their simultaneous eigenmatrices. In particular, in the presence of a continuous weak symmetry in Eq. (3), there exists a weakly symmetric representation such that withH andJ j supported as in Eq. (9), but with e i(φ k −φ l ) = e iδ j replaced by s k − s l = δ j , where s k is an eigenvalue of S corresponding to an eigenspace H k .

B. Construction
We now give two constructions of weakly symmetric representations from a given representation of the master operator in Eq. (1), that is, a Hamiltonian H and a set of jump operators {J j } j . In the first construction, the Hamiltonian is projected on the symmetric eigenspace of the symmetry superoperator, while the jump operators are projected on all its eigenspaces, so that their number n in general increases to m-fold, where m is the number of distinct eigenvalues of the symmetry superoperator. Here the knowledge of symmetry eigenspaces is assumed, which generally requires diagonalizing the dim(H) × dim(H) matrix of the symmetry operator U (see Fig. 2). In the second construction, the number of jump operators does not increase as a weakly symmetric representation with the minimal number of jump operators is constructed, but at the cost of diagonalizing two matrices of size n × n (see Fig. 3).

Weakly symmetric representation by dynamical decoupling
In this construction we use the fact that for dynamics with a weak symmetry in Eq. (2), The right-hand-side limits in Eqs. (11) corresponds to the projection of the master operator on the symmetric part under its transformations U (·)U † , and in the construction we consider this projection applied to individual terms appearing in the master equation in Eq. (1) (see Fig. 2). Note that such a projection occurs when the dynamics is in fact composed of the system master dynamics governed by L and much faster unitary dynamics corresponding to U , in which case the former acts as a perturbation of the latter (see Supplemental Material of Refs. [26,27]). Therefore, a weak symmetry can be facilitated by dynamical decoupling [28,29] at a rate much  (19)] gives weakly symmetric jump operators [Eqs. (20)].
faster than system dynamics but much slower than relaxation of the environment (see Refs. [30,31]).
Step 1. A symmetry operator U is diagonalized to find its eigenspaces.
Step 2. Eigenspaces of the symmetry superoperator U are constructed.
Step 3. The Hamiltonian and jump operators are replaced by their projections on U eigenspaces [cf. Eq. (9)]: with e iφ k denoting an eigenvalue of U corresponding to an eigenspace H k and e iδ denoting an eigenvalue of U . Construction for Abelian weak symmetries. Considering Eq. (11) for all Abelian weak symmetries present, Eq. (12a) holds with subspaces H k defined as intersections of eigenspaces of the symmetry operators, while the sum in Eq. (12b) is replaced k, l with the ratios of the symmetry operators eigenvalues corresponding to a given set of eigenvalues for the corresponding symmetry superoperators. In particular, for a continuous symmetry in Eq. (3), and an eigenspace H k of S corresponding to an eigenvalue s k , the Hamiltonian is constructed as in Eq. (12a), while a jump operator J j is replaced by the set Proof of Eq. (12). First, note that U L U † corresponds to the master operator with the Hamiltonian U (H ) and the jump operators U (J j ), as seen, for example, in the Liouville representation. The master operator L in Eq. (5) is linear in H ⊗ 1 and J † j J j ⊗ 1, so that in the limit of the right-hand side in Eq. (11) they are replaced by their projection on the symmetric subspace of U ⊗ I, where I is an identity superoperator, that is,H ⊗ 1 of Eq. (12a) and j,e iδ ⊗ 1. Similarly, the term J j ⊗ J * j is projected on the symmetric subspace of U ⊗ U * and is thus replaced by e iδJ j,e iδ ⊗J * j,e iδ .

Minimal weakly symmetric representation
The steps of this construction are motivated by following two facts (see Fig. 3). First, a representation of the master operator with the traceless Hamiltonian and orthogonal traceless jump operators is uniquely defined (up to degeneracy in jump rates) and corresponds to the minimal number, n min n, of jump operators (see Ref. [23]). Second, the set U (H ), {U (J j )} j is a representation of UL U † [32], and thus, in the presence of the weak symmetry, it is also a representation of L [7]. Since U does not change the orthogonality and trace of the operator being transformed, we have that in the presence of weak symmetry, the traceless Hamiltonian is necessarily symmetric, while orthogonal traceless jump operators with the same rate are transformed unitarily by U , and thus they can be chosen as its eigenmatrices.
Step 1. Traceless jump operators are constructed by introducing while the Hamiltonian is replaced bỹ in order to leave the master operator in Eq. (1) unchanged [a further shift of the Hamiltonian by −Tr(H )1 only introduces a global phase, and is also symmetric, and thus will be omitted].
Step 2. A Hermitian matrix of the scalar products, 2 is diagonalized in order to define via its orthonormal eigenvectors, with the rate determined by the corresponding eigenvalue We reorder jump operators with decreasing λ j and neglect j with λ j = 0, as J j = 0 ( j n min ).
Step 3. For a weak symmetry in Eq. (2), the Hamiltonian in Eq. (14) is symmetric, while the set of orthogonal jump operators in Eq. (16) is is a unitary matrix 3 which is block diagonal in the eigenspaces of . The jump operators determined by the orthonormal eigenvectors of U, U u j = e iδ j u j , are eigenmatrices of the symmetry superoperator, j = 1, ..., n min . In particular, when diagonalizing blocks in U, the jump operators in Eq. (20) are chosen orthogonal, as 3 U is unitary as U † corresponds to the superoperator U † , which generates the unitary transformation U φ = e iφS of orthogonal jump operators under U φ = e iφS [cf. Eq. (19)] and is block diagonal in the eigenspaces of [Eq. (17)].

C. Nonuniqueness
In Sec. III B we constructed two generally different representations of the master equation, which proves that a weakly symmetric representation [Eq. (8)] is nonunique. Here, we characterize the freedom in the choice of weakly symmetric representations.
A general weakly symmetric representation with traceless jump operators is described by an n × n min isometry V [(V † V) jk = δ jk ] that does not mix U eigenspaces [(V) jk = 0 and (V) jl = 0 take place only for e iδ k = e iδ l ], that is, the set of jump operators {J j } j defined in Eq. (20) is replaced by { n min k=1 (V) jkJk } n j=1 (cf. Ref. [23]). In that case, jump operators are generally not orthogonal, with the scalar product between the jth and kth jump given by ( Eqs. (15) and (17)]. A general weakly symmetric representation features to shifted symmetric jump operators (cf. Ref. [23]). That is, symmetric jump operators in a weakly symmetric representation need not to be traceless, as can be shifted by a constant since [U, 1] = 0. Indeed, any jump operatorJ j with e iδ j = 1 in a weakly symmetric representation can be replaced byJ j + a j 1 with a j ∈ C, while the symmetric HamiltonianH is transformed toH + b1 − i j: e iδ j =1 (a * jJ j − a jJ † j )/2 with b ∈ R [cf. Eqs. (13) and (14)].

IV. QUANTUM TRAJECTORIES WITH WEAKLY SYMMETRIC REPRESENTATIONS
We now briefly discuss implications of the presence of a weak symmetry for the structure of quantum trajectories and the survival of coherences between symmetry eigenspaces. This structure simplifies the construction of the master operator and reduces both the memory and processing required for QJMC simulations, as we explain in Sec. V.

B. Simplified quantum trajectories
We now discuss how quantum trajectories simplify for a weakly symmetric representation (see Fig. 4).

Symmetric initial states
Consider the system initially supported in a symmetry eigenspace H l , U |ψ 0 = e iφ l |ψ 0 , which we refer to as a symmetric initial state since U (|ψ 0 ψ 0 |) = 0. For a weakly symmetric representation, the system state remains supported in the same eigenspace when no jumps take place [cf. Eq. (24b) and see Fig. 4

(a)],
because of the symmetry of the effective Hamiltonian, Occurrence of the first jumpJ j 1 at time t 1 transforms the state |ψ t 1 from H l into the symmetry eigenspace H k with the eigenvalue e iφ k = e i(δ j 1 +φ l ) [cf. Eqs. (9b) and (24a)], Only for a symmetric jump, e iδ j 1 = 1, the symmetry eigenspace remains unchanged, H k = H l . The system state remains in H k until the next asymmetric jump [see Fig. 4(a)]. Therefore, for an initially symmetric system state, quantum trajectories with a weakly symmetric representation are In contrast, in a general representation, coherences between symmetry eigenspaces are present in individual quantum trajectories but interfere destructively in the average of Eq. (23). This illustrates the fact that quantum stochastic dynamics with a weakly symmetric representation features the weak symmetry, as do the dynamics of a density matrix under the effective Hamiltonian or the action of any of jump operators [cf. Eqs. (28) and (30)]. The weak symmetry is then inherited by the average dynamics [see Fig. 4(b)]. Analogous results hold for the case of Abelian weak symmetries.

General initial states
For an initial state in a superposition of symmetry eigenstates, |ψ 0 = l 1 H l |ψ 0 , the coherences between symmetry eigenspaces are in general maintained in a quantum trajectory, since in a weakly symmetric representation asymmetric jump operators connect many pairs of symmetry eigenspaces [see Eq. (9b)]. Nevertheless, no new coherences with respect to the eigenvalues of the symmetry superoperator are created, since the action of individual jump operators also features the weak symmetry [cf. Fig. 1(b)]. Furthermore, for a unique stationary state, from the ergodicity [Eq. (27)], any coherences in a quantum trajectory must either interfere destructively in the time average or decay to 0 over time. When not only the asymptotic time average in a quantum trajectory but also the asymptotic time distribution is independent from an initial state, as is the case for generic dynamics [34], coherences in quantum trajectories necessarily decay to 0.

V. NUMERICAL APPLICATIONS OF WEAKLY SYMMETRIC REPRESENTATIONS
We now explain how weakly symmetric representations can be used to simplify a construction of the master operator as well as QJMC simulations.

A. Simplified master operator construction
As discussed in Sec. II B 2, the master operator with a weak symmetry becomes block diagonal in the Liouville representation when using a basis corresponding to symmetry eigenspaces. This simplifies both its diagonalization and numerical integration. We now show how its construction in such a basis can be further simplified by using a weakly symmetric representation.
We have [cf. Eqs. (5) and (25)] For dynamics with a weak symmetry in Eq. (2), we consider a weakly symmetric representation in Eq. (8). From Eq. (9), a subspace H k ⊗ H l is mapped onto itself only by the symmetric effective Hamiltonian and symmetric jump operators (cf. Fig. 5): with Fig. 5): The construction of the master operator in the Liouville representation with a weakly symmetric representation is made efficient in two ways. First, only the nontrivial action of L within the eigenspaces of the symmetry superoperator U  33) and (34). For a discrete weak symmetry with U with trivially degenerate eigenvalue ratios (or a continuous weak symmetry with S with trivially degenerate gaps), H k ⊗ H l with k = l is always mapped onto itself [dotted gray block vanishes; cf. Fig. 1(b)].
is computed, i.e., only nonzero blocks of L are constructed. 5 Second, transformations between individual subspaces H k ⊗ H l can be computed using only the operators in the representation that correspond to a specific eigenvalue of U ; see Fig. 5. Therefore, for a minimal symmetric representation, the maximal number of jump operators that contributes in Eq. (33) is the maximal number of symmetric jump operators, k dim(H k ) 2 − 1, while in Eq. (34) it is the maximal number of jump operators with the eigenvalue e i(φ k −φ k ) , i.e., l,l : e i(φ l −φ l ) =e i(φ k −φ k ) dim(H l ) dim(H l ) (as can be seen, for example, by using the Choi representation [23]).
The simplified construction of the master operator is a direct consequence of the simplified structure of quantum trajectories discussed in Sec. IV. Indeed, an infinitesimal change in the average system state is a consequence of a change in quantum trajectories, either due to the symmetric effective Hamiltonian [cf. Eq. (28)] or to an occurrence of a jumpJ j [Eq. (30)]. The mapping of H k ⊗ H l for k = l is then determined by quantum trajectories originating in H k , while for k = l it is determined by quantum trajectories of superpositions in H k ⊕ H l (cf. Fig. 5).
Analogous results hold for Abelian weak symmetries. In particular, for a continuous weak symmetry in Eq. (3) and a weakly symmetric representation in Eq. (10), we have Eq. (33) with e iδ j = 1 replaced by δ j = 0, and Eq. (34) with e iδ j = e i(φ k −φ k ) replaced by δ j = s k − s k for H k ⊗ H l is mapped onto H k ⊗ H l within the same eigenspace of S (that is, for s k − s k = s l − s l ).

B. Simplified QJMC simulations
The QJMC approach (see, e.g., Refs. [5][6][7][8][9]) is used to generate quantum trajectories in order to obtain dynamics of the average system state via the empirical mean [cf. Eq. (23)] or stationary states of the system [by considering trajectories longer than the relaxation time, or, in the case of a unique stationary state, via Eq. (27)]. The advantage of this method in comparison with solving Eq. (1) lies in considering linear operators on pure states in the system space H rather than the master operator on density matrices in the space isomorphic to H ⊗ H (cf. Fig. 1). Since, by definition, the average dynamics governed by the master operator in Eq. (1) does not depend on its representation, it can be chosen to simplify QJMC simulations. Here, we discuss how for the dynamics with a weak symmetry, weakly symmetric representations can be utilized (see also Refs. [7,9,35]).

Algorithm
In order to construct a quantum trajectory up to a finite time t, each step consists of two parts. First, for an initial state |ψ 0 , time t 1 of the first jump is found by drawing a random uniformly distributed number u 1 ∈ [0, 1], which represents the probability of no jump occurring until t 1 [cf. Eq. (26b)], Second, if t 1 > t, the normalized quantum trajectory at time t is given by normalized Eq. (24b). Otherwise, a jump takes place of the type j 1 drawn with the probability p j 1 proportional to its instant rate, and the system state is updated as |ψ t 1 +dt 1 ( j 1 , t 1 ) = J j 1 |ψ t 1 [cf. Eq. (24a)]. Then, in order to find time t 2 of the next jump, the step is repeated with |ψ 0 replaced by the normalized conditional state |ψ t 1 +dt 1 ( j 1 , t 1 ) / |ψ t 1 +dt 1 ( j 1 , t 1 ) and t 1 replaced by time t 2 − t 1 between the first and the second jumps. This is done until t n+1 > t, in which case the normalized quantum trajectory after the occurrence of n jumps is given at time t by normalized Eq. (24a). The main computational difficulty is the evaluation of jump occurrence times [Eq. (35)], which requires finding a norm of a system state evolving under the non-Hermitian effective Hamiltonian H eff . Therefore efficient computation of the no-jump dynamics is necessary, e.g., via the exact diagonalization of H eff or via fourth-order Runge-Kutta integration [8,9]. We note that, alternatively, a discretization of quantum trajectories to finite but small time steps δt can be considered. Here, at each step, occurrence of a jump is decided when a (c) This state is normalized and used as the initial state in the next step, which is repeated until t n+1 > t when the state at time t is chosen as normalized Eq. (24a), inside the symmetry eigenspace with the eigenvalue e iδ jn · · · e iδ j 1 e iφ l . random uniformly distributed number [0,1] is smaller than δt j ψ|J † j J j |ψ , which approximates, up to the linear order, the probability of a jump occurring for the system in |ψ at the beginning of the step [cf. Eq. (35)], upon which the type of jump j is drawn according to its rate [cf. Eq. (36)]. The system state is then updated to J j |ψ or, when no jump occurs, to (1 − iδtH eff )|ψ , and subsequently normalized before the next step. This originally proposed approach [6,7] is equivalent to first-order Euler integration of Eq. (35) [8]. Nevertheless, another important factor remains: manipulating the conditional system state generally described by dim(H) coefficients [for example, in Eq. (36) or when updating and normalizing the state].

Simplified algorithm
We now explain how the complexity of the QJMC algorithm can be lowered thanks to a weak symmetry by considering a weakly symmetric representation. Simulations are simplified, in particular, for symmetric initial states, which is relevant for system dynamics with a unique stationary state and for averages of symmetric system observables (see Fig. 6).
Symmetric initial states. For a symmetric initial state, any quantum trajectory is confined to only a single symmetry subspace at a time [cf. Eq. (31) and see Fig. 4(a)]. Therefore the system state is described at any time by at most max k dim(H k ) coefficients (and the label k for the occupied eigenspace H k ). Furthermore, the evolution with the symmetric effective Hamiltonian in each step of the algorithm can be integrated, if needed, solely on the currently occupied symmetry subspace [that is, for H l , 1 H lH eff 1 H l in Eq. (35)]. Furthermore, the stochastic dynamics remains exactly the same upon replacing each of the jump operators J j with the set {1 H kJ j 1 H l } k,l: e iδ j =e i(φ k −φ l ) [cf. Eq. (30) and Fig. 4], with the rates of jumps for a state in H l simplified While these improvements are known for the open quantum dynamics microscopically defined by a weakly symmetric representation corresponding to a continuous weak symmetry [9], such as particle losses in a many-body system corresponding to the U (1) symmetry generated by the particle number [35], in this work we show how to implement them for any weak symmetry by constructing a weakly symmetric representation.
General initial states. For a general initial state, |ψ 0 = l 1 H l |ψ 0 , in each step of the algorithm, the effective Hamiltonian can be integrated, if needed, independently in each symmetry subspace, (35)]. Furthermore, a type of occurring jump can be determined with respect to the sum of its rates in individual subspaces, Eq. (36)], and the updated state corresponds to the superposition of the updated states, We conclude that QJMC simulations with weakly symmetric representations are simplified in a way comparable to the strong symmetry case with a Hamiltonian and all jump operators being symmetric [16,17]. Indeed, in that case there exists a stationary state ρ (k) ss inside each symmetry eigenspace H k , which can be obtained from quantum trajectories for an initial state within that subspace evolving with the effective Hamiltonian 1 H k H1 H k and jump operators {1 H k J j 1 H k } j . Similarly, for a general initial state, the dynamics can be solved independently in each symmetry eigenspace.

C. Sparsity
Hamiltonians and jump operators for many-body system involving only few-body terms are sparse in any basis composed as a tensor product of local bases. Since the number of few-body operators scales linearly in the system size, rather than exponentially as the dimension of the system Hilbert space, the master operator is also sparse [cf. Eq. (5)]. This allows for a significant computational speedup in its diagonalization or numerical integration.
Although a weakly symmetric representation features a Hamiltonian and jump operators which are linear combinations of the original operators (cf. Secs. III B and III C), their number scales linearly with system size and thus does not significantly affect the sparsity. In order to exploit the weak symmetry, either for construction and diagonalization of the master operator or for QJMC simulations, however, it is necessary to work in the basis of eigenspaces of the symmetry operator U (cf. Figs. 1, 5, and 4). For a local symmetry (U being a tensor product of local unitaries), e.g., a global rotation of spin systems, the basis of symmetry eigenspaces can be chosen separable. For a nonlocal symmetry, such as a translation symmetry, symmetric states are entangled. Therefore, similarly as for symmetries in closed quantum dynamics, it needs to be judged on a case-by-case basis whether exploiting a weak symmetry leads to improved numerical simulations of the open quantum dynamics (see Fig. 7 and Sec. VI).

VI. EXAMPLES
Finally, we give closed formulas for weakly symmetric representations of many-body open quantum system dynamics with translation and rotation symmetries. We utilize these representations for numerical simulations of a spin-1 chain with nearest-neighbor interactions and local dissipation shown in Figs. 7 and 8, which features both symmetries.
where q denotes a quasimomentum and T (J q,α ) = e i 2πq A symmetric basis can be determined as follows. For a local basis of the system chosen as product states of identical subsystem bases, each element belongs to a cycle of a length l under the translation symmetry, where l divides N. The planewave superpositions of the basis elements with quasimomenta q l corresponding to that cycle are eigenstates of T corresponding to eigenvalues e i2πq l /l , where q l = 0, 1, ..., l − 1, so that the effective quasimomenta is q = Nq l /l. Therefore the dimension of a symmetry subspace indexed by q equals the number of cycles with lengths l = N/d q , where d q is a common divisor of q and N [cf. Similarly, collective plane-wave jump operators in Eq. (37) being linear combinations of N jump operators feature at most N 3 z α nonzero elements in each column, where z α is the maximal number of nonzero entries in columns of the matrix in the local basis for a jump of type α. For a local Hamiltonian H and local jumps J j,α , z and z α are independent from N. these generators do not commute with S z . In contrast, when the weak translation symmetry is also present, e.g., λ j ≡ λ, V jk ≡ V | j−k| , and W jklm ≡ W | j−k||l−m| (see Fig. 8), the jump operators in Eq. (37) for α = z, +, −, are simultaneously eigenmatrices of S z and T (note that [S z , T ] = 0). Note that the jump operatorsJ q,± with q < N are uniquely determined (cf. Sec. III C).

VII. CONCLUSIONS
In this article, we investigated how Abelian weak symmetries in Markovian dynamics of open quantum systems can be translated into their quantum stochastic dynamics. We showed how to construct weakly symmetric representations of a master operator governing the dynamics, for which quantum trajectories are symmetric, i.e., are found within a single symmetry eigenspace at a time, whenever initial system states are chosen symmetric. This enabled us to exploit weak symmetries of the dynamics for simplifying the QJMC algorithm, with the memory and processing required for simulations reduced in a way akin to the case of strong symmetries. We also demonstrated how the efficiency in constructing the Liouville representation of the master operator can be improved-a result directly relevant for solving the system evolution via diagonalization or numerical integration of the master operator. Finally, we note that the existence of weakly symmetric representations does not rely on the dynamics being time independent, as considered, for simplicity, in this work. For a time-dependent dynamics that features the same weak symmetry at all times, a time-dependent weakly symmetric representation can be analogously constructed from given time-dependent Hamiltonian and jump operators.