Hilbert Space Fragmentation in Open Quantum Systems

We investigate the phenomenon of Hilbert space fragmentation (HSF) in open quantum systems and find that it can stabilize highly entangled steady states. For concreteness, we consider the Temperley-Lieb model, which exhibits quantum HSF in an entangled basis, and investigate the Lindblad dynamics under two different couplings. First, we couple the system to a dephasing bath that reduces quantum fragmentation to a classical one with the resulting stationary state being separable. We observe that despite vanishing quantum correlations, classical correlations develop due to fluctuations of the remaining conserved quantities, which we show can be captured by a classical stochastic circuit evolution. Second, we use a coupling that preserves the quantum fragmentation structure. We derive a general expression for the steady state, which has a strong coherent memory of the initial state due to the extensive number of non-commuting conserved quantities. We show that it is highly entangled as quantified by the logarithmic negativity.

FIG. 1. Schematic representation of the setup.The Hilbert spaces of both the pair-flip (PF) and Temperley-Lieb (TL) models fragment into exponentially many Krylov subspaces (solid filled blue and green squares respectively).The degenerate Krylov subspaces of the TL model are contained in the same grey squares.The dephasing noise Lj = S z j connects some of the fragmented subspaces of the TL model (onsite dissipative coupling in blue), such that the fragmentation reduces to the classical one of the PF model.Nonetheless, the quantum fragmentation is preserved when using specific two-site dissipative couplings.
First examples of fragmentation in an entangled basiswhich we denote quantum fragmentation (QF)-have been only recently proposed [46].Reference [46] put forward an algebraic approach using the mathematical notion of bond and commutant algebras to characterize the set of conserved quantities, which also provides a systematic way to explore the differences among these two types of fragmentation.
In realistic settings, quantum many-body systems are never perfectly isolated from their surrounding environment.This raises the question to which extent the phenomena related to HSF -in particular QF, which takes place in an entangled basis-are affected by couplings to a bath.Note that in case of MBL, the localization is destroyed when the system is locally coupled to a dissipative bath [47][48][49][50].In Ref. [51], CF in open systems due to weak symmetries was studied by exploiting the resulting integrable structure, allowing to obtain the spectrum of the Lindbladian within all invariant subspaces.Here we aim to provide an understanding of the generic behavior that quantum fragmented models can display in the presence of a dissipative bath and consider a family of QF models introduced in Ref. [46].Building up on previous works on the stationary state structure for Lindbladian evolution in the presence of (conventional) conserved quantities [52][53][54], we investigate systems described by the commutant algebra formalism [46] and focus on different strong symmetries [53,54], i.e., symmetries preserved by both the Hamiltonian and every jump operator.We start by considering a dephasing coupling [55], where the system locally couples to a bath that we find to eventually reduce the QF to CF.Still, the system preserves a large amount of information of the initial state even at infinite time due to the extensive degeneracy of stationary states.On the other hand, a dissipative environment can be engineered and exploited to create exotic non-equilibrium dynamics [56][57][58][59][60].For example, Ref. [58] proposed to efficiently drive the system to the desired pure state as the unique stationary state by engineering dissipative couplings.Motivated by this, we consider a fine-tuned coupling which preserves the structure of the QF system.Interestingly, we find that the system evolves to a highly entangled stationary state.Moreover, we propose this as a simple protocol to decide whether a system is quantum fragmented.
The remainder of the paper is organized as follows.In Sec.II, we briefly review the commutant and bond algebras formulation for isolated fragmented systems [46], and generalize it to open quantum systems focusing on strong symmetries.In Sec.III, we then introduce two related fragmented models, the Pair-flip (PF) and the Temperley-Lieb (TL) model, which exhibit classical and quantum fragmentation, respectively.We study the TL model under dephasing noise in Sec.IV, which leads to a breakdown of quantum fragmentation to the classical fragmentation of the PF model.In Sec.V, we couple the TL model to the structure-preserving noise, which preserves the original quantum fragmentation of the TL model.We analytically derive the stationary states of the dynamics under both couplings that we use to predict saturation values of two-point correlators and two different entanglement measures, the logarithmic negativity, and the operator space entanglement, and compare them with numerical simulations.We conclude in Sec.VI by summarizing our main findings and discussing open questions.Finally, we consign more technical aspects of our work to the appendices.
The phenomenon of HSF arises as a consequence of certain constraints being imposed on the dynamics of many-body systems.Given a family of Hamiltonians H = j J j h j parameterized by real coefficients {J j }, fragmentation is a property that is completely characterized by the local terms {h j } and thus holds for any choice of coefficients.This distinguishes HSF from other symmetries that might appear for certain choices of J j , such as translation invariance with uniform J j .Reference [46] formalized this observation using the language of bond and commutant algebras for isolated quantum systems, which we will review in the following.A bond algebra A is the algebra generated by arbitrary linear combinations of products of the local terms {h j }, together with the identity operator 1.The corresponding commutant algebra C is the set of conserved quantities, namely the centralizer of A including all operators that commute with every local term We refer to the latter using the shorthand notation C = {h j } .Both A and C are von Neumann algebras, i.e., they include the identity operator and are closed under conjugation [61].Importantly, every element in C commutes with every element in A, i.e., they are the centralizers of each other.As such the Hilbert space can be decomposed into irreducible representations of C × A [61,62], where λ and H (A) λ are the d λ and D λ dimensional irreducible representations of C and A, respectively.This decomposition implies that the elements of the bond algebra h A ∈ A, generate independent dynamics within H (A) λ while acting trivially on H (C) λ .Therefore, for fixed λ, there are d λ degenerate Krylov subspaces or fragments with dimension D λ .We will denote the degenerate Krylov subspaces as K λ α , with α = 1, ..., d λ , and omit λ if there is no degeneracy.
The formulation in terms of bond and commutant algebras provides a unifying framework to describe the decomposition of the Hilbert space, which applies to both conventional and unconventional symmetries like HSF [46,63].The difference appears in the scaling of the dimension of the commutant dim(C) = λ d 2 λ with system size: It scales exponentially for fragmented systems while at most polynomially for conventional symmetries.When the commutant C is Abelian, every irreducible representation is one-dimensional (d λ = 1) and hence, the Hilbert space reduces to a direct sum of non-degenerate Krylov subspaces, H = α K α .Projectors Π α = β |ψ αβ ψ αβ | onto those subspaces span the commutant, where {|ψ αβ } is an orthonormal basis in K α .On the other hand, non-Abelian commutants include larger dimensional irreducible representations d λ > 1, corresponding to degenerate Krylov subspaces.In this case, the projectors Π λ α onto different Krylov subspaces K λ α span a maximal Abelian subalgebra of C, while the full C is generated by not only the projectors but also the intertwine operators between degenerate ones, Π λ αα = β |ψ λ αβ ψ λ α β | [46].For example, the commutant algebra of SU(2)-symmetric systems is non-Abelian and contains non-commuting conserved quantities such as S x tot , S y tot , and S z tot .The total spin representation λ is given by the eigenvalues of ( S tot ) 2 as λ(λ + 1).There are d λ = 2λ + 1 degenerate Krylov subspaces with the same λ, which are labeled by different spin-z projections S z tot = −λ, −λ + 1, ...λ, leading to the Hilbert space decomposition as in Eq. ( 2) [64].
Fragmentation can be classified as either classical or quantum.A system is said to be classically fragmented if one can find a common eigenbasis of product states for all elements in a maximal Abelian subalgebra of the commutant.This means that the Krylov subspaces can be spanned by a product state basis.Otherwise a system is said to be quantum fragmented.By this definition, we associate CF with the existence of a basis of product states and QF with an entangled basis-which is different from the commutant being Abelian or not.Specifically, an entangled basis can also appear for Abelian commutants.For example, for an SU(2)-symmetric system, adding the term S z tot preserves the Hilbert space structure in an entangled basis but breaks the degeneracy of the Krylov subspaces, leading to the so-called dynamical SU(2) symmetry and an Abelian commutant [63].Note that the current definition of QF is still not ideal since it leaves some room for ambiguous or trivial examples.

B. Lindblad dynamics of fragmented systems
We study the dynamics of fragmented systems coupled to a Markovian bath described by a Lindblad master equation, dρ dt = L(ρ) (see Fig. 1).Here L is a Liouvillian superoperator with [55,65] where the positive coefficients γ j correspond to the decay rates, {L j } are jump operators describing the coupling to a bath, and we set = 1.Equivalently, the time evolution of an operator in the Heisenberg picture is generated by the adjoint of the Liouvillian superoperator, dO dt = L † (O).Of particular interest to us is the stationary state as an eigenstate of L with zero eigenvalue, namely L(ρ ss ) = 0. Generally, for a Liouvillian without symmetries, there is a unique stationary state and it preserves no information of the initial state.On the other hand, symmetries and conserved quantities can lead to multiple stationary states and a memory effect in the long time limit [52][53][54]66].A simple case is the presence of a strong unitary symmetry S that is preserved by both the Hamiltonian and every jump operator, i.e., [S, H] = [S, L j ] = [S, L † j ] = 0 for all j [53,54].The space of bounded operators B(H) decomposes into orthonormal subspaces, B αα = span{|ψ α ψ α |}, where |ψ α is an eigenstate of S with eigenvalue s α .Each subspace labeled by different quantum numbers of S evolves independently since LB αα ⊆ B αα .Thus, the stationary state inherits the block diagonal structure given by the symmetry, which leads to at least as many distinct stationary states as the number of symmetry sectors [53,54].
Let us now investigate the phenomenon of HSF in the strong symmetry sense.In Lindblad systems, the dynamics is generated by the Hamiltonian H = j J j h j and the jump operators {L j }.Reference [67] considered the commutant H, {L j }, {L † j } associated with the (total) Hamiltonian and the jump operators, which was shown to give a complete set of conserved projectors onto mutually orthogonal subspaces with independent dynamics in B(H).Note however that the analysis only applies when the conserved operators form an algebra1 .In the language of Ref. [63], H, {L j }, {L † j } corresponds to a local algebra rather than to a bond algebra, since H is an extensive sum of local terms.To extend the analysis of HSF in terms of bond and commutant algebras to open quantum systems, we define the open bond algebra A O = {h j }, {L j } where we focus on Hermitian L j , and the corresponding open commutant = 0 for all j, such that each subspace evolves independently.Note that these Π λ α project onto minimal (irreducible) subspaces of the dynamics generated by L, as they span the maximal Abelian subalgebra of the open commutant [67].All together we find that the operator space B(H) decomposes into orthogonal, invariant, minimal subspaces B αα [67,68].As stated above, the existence of non-unique stationary states is now guaranteed by these subspaces, where now the degeneracy of the stationary state scales exponentially with system size due to HSF.

III. MODEL AND SETUP
We study the dynamics of quantum fragmented systems coupled to a dissipative environment by considering the family of Temperley-Lieb (TL) models as a concrete example.
First, we introduce the closely-related spin-1 pair-flip (PF) model, which exhibits CF, i.e., it is fragmented in a product-state basis.The Hamiltonian is given by where α, β denote different spin-z components {−, 0, +}, and g αβ j,j+1 and l jα are arbitrary real coefficients.We assume open boundary conditions (OBC) and even number of sites for convenience.The constrained dynamics of the PF model can be visualized by mapping product states in the computational basis to colored pairs and dots.Specifically The PF model has two independent U(1) charges, which are given by N + = j (−1) j N + j and N − = j (−1) j N − j , with N α j = (|α α|) j .These U(1) symmetry sectors further split into smaller Krylov subspaces labeled by a nonlocal invariant [69], which we will discuss in the following.Starting from a product state, we first connect all the adjacent spins with the same color from left to right.Next we remove the paired spins and repeat the first step until there are only unpaired spins with a different color from their nearest neighbors to the left.The unpaired spins are then referred to as dots.We denote dot patterns of size 2λ as A λ .Let us for example consider a state with the dot pattern ( ), We observe that dot patterns, i.e., the color and sequence of unpaired spins, are invariant under the action of a pair-flip, providing non-local and mutually commuting conserved quantities similarly to Ref. [19].Thus, each Krylov subspace can be labeled by a dot pattern.
The number of different dot patterns grows exponentially with system size and thus the Hilbert space fragments into exponentially many Krylov subspaces in the local z basis.Since the fragmentation occurs in a basis of product states -common eigenbasis of all elements of the Abelian commutant-the PF model exhibits CF.See schematic representation appearing in Fig. 1.
As a special case of the PF model, the TL model is at least as (classically) fragmented as the PF model.In addition to the SU(3) symmetry, the constrained dynamics conserves extended dot patterns, including colored dot patterns of the PF model, e.g., | j,k , and additional entangled dot patterns, e.g., [46].Note that the choice of the dot states is not unique due to the non-Abelian nature of C. The following example shows that the dot pattern | j,k is conserved: The TL model is then block-diagonal in an entangled basis given by the dimers and dots configurations (more examples are shown in App.A 1). Thus the TL model exhibits QF.Moreover, the resulting commutant algebra C TL is non-Abelian and the dimension of the irreducible subspaces are d λ ≥ 1 [46,71].Therefore, there are d λ degenerate Krylov subspaces for fixed λ, which are labeled by different dot patterns with the same 2λ length.Note that in the previous discussion, we distinguished between CF and Abelian commutant, as well as between QF and non-Abelian commutant.For example, one can find systems with an Abelian commutant which nonetheless require an entangled basis [63].
Following the distinction between strong and weak fragmentation as discussed in Ref. [16], we verify that both the PF model and the TL model exhibit strong fragmentation with respect to the full Hilbert space.The dimension of the largest Krylov subspace scale as D max /3 N ∼ exp(−aN ) with a < 1. See App.A 2 for additional details.
To study the effect of fragmentation on the Lindblad evolution, we discretize the dynamics and implement a local random quantum circuit including both Hamiltonian Liouvillian gate Unitary gate FIG. 2. Lindblad random circuits.A single time step for random circuit evolution for (a) a closed system with twosite unitary gates Uτ,j = e −iJ j h j,j+1 ; and (b) the Lindblad evolution with Liouvillian gates Uτ,j = e L j,j+1 , where {Jj} are random coefficients extracted from a uniform distribution Jj ∈ [0.8, 1.2].Here, the blue circles represent the initial state or density matrix.
and Lindblad evolutions.This implementation breaks energy conservation and translation symmetry, and only preserves those quantities belonging to the commutant algebra.The setting of random circuits is shown in Fig. 2. Every time step includes two consecutive layers of nonoverlaping gates, such that after t time steps the time evolution is given by The Liouvillian gates are superoperators given by U τ,j = e Lj,j+1 (see Fig. 2b), with where {J j } are uniformly distributed random coefficients for different sites j and time steps τ .The dissipation term is ) for one-site jump operators, and D j,j+1 (ρ) = γ(L j,j+1 ρL † j,j+1 − 1 2 {L † j,j+1 L j,j+1 , ρ}) for two-site jump operators.In the following, we use the dimensionless hoppings J j to be uniformly distributed in the interval [0.8, 1.2].Moreover, we choose γ j = γ as this does not affect our results [68].When γ = 0, the Liouvillian gates become random unitaries with the overall phase fluctuating around π.This implementation allows us to compare our numerical results with the analytic prediction obtained using the formalism introduced in the previous section.

IV. DEPHASING NOISE
We first consider a dephasing noise given by L j = S z j .For many-body localized systems, such coupling delocalizes the system and drives it to an infinite temperature state ρ ∝ 1 [47][48][49][50].For the TL model, however, the dephasing noise preserves the CF while breaking the QF.When considering the whole Hilbert space, this turns into non-ergodic behavior and extensively degenerate stationary states.
The mechanism for the breakdown from the QF of the TL model to the CF of the PF model is shown in Fig. 1.Intuitively, the TL model is symmetric with respect to different color pairs due to the SU(3) symmetry, while the dephasing noise distinguishes different colors.However, this respects the CF, as the jump operators are elements of the PF bond algebra, S z j ∈ A PF .Moreover, any element of this algebra can be written as linear combinations of products of elements in the TL bond algebra and the dephasing jump operators as explicitly shown in App.B 1. Therefore, the corresponding bond algebra is given by the PF one {h j }, {S z j } = A PF with open commutant C O = C PF .This implies that the symmetries of the Liouvillian are those of the PF model.In this section, we study the effect of the breakdown of quantum fragmentation, and sketch the derivation of the stationary state for this case.

A. Stationary states with classical fragmentation
We now derive the stationary state of the TL model under dephasing noise.As we just showed, both the Hamiltonian and the jump operators preserve the CF of the PF model, with [h j , Π α ] = [L j , Π α ] = 0, ∀j.Therefore, the operator space is decomposed into orthogonal subspaces with independent dynamics, or equivalently, LB αα ⊆ B αα , where we denote the diagonal subspaces B α ≡ B αα .This is the natural extension of strong symmetry for fragmented systems.
Next, we show that there is a unique stationary state within each B α .As all jump operators are Hermitian, the infinite temperature state ρ ∝ 1 is a stationary state in B(H).Therefore, there exists a stationary state This is because the dissipation induces full decoherence within each invariant subspace.Moreover, as the projectors {Π α } span the maximal Abelian subalgebra of the open commutant [46], these invariant subspaces B α are minimal subspaces [67].Therefore, the stationary state within each B α is unique [67].Additional details can be found in App.B 2.
Combining the stationary state structure (i.e., a unique stationary state within each minimal subspace) and the corresponding conserved quantities {Π α }, we find that the general expression of the stationary state is given by The coefficients c α ∈ R are the weights of the initial state within the diagonal subspaces K α .The stationary state preserves the weight c α , while all the off-diagonal (coherent) information is lost.A more detailed derivation can be found in App.initial states scales exponentially with the system size, signalling a strong memory effect.In the following, using the expression in Eq. ( 13), we analyze the long-time behavior of the TL model under dephasing noise.Figure 3 shows an example of the evolution of the density operator (written in the local z-basis) by exact diagonalization (ED).To compare with the case of quantum fragmentation we use the initial state which has non-zero overlap only with three Krylov subspaces: the fully-paired subspace (with zero dots) and other two labeled by the dot patterns ( ) and ( ).At long times, all off-diagonal matrix elements vanish.The stationary state is then the direct sum of projected identities within the diagonal blocks, with the weight determined by the initial state.

B. Infinite temperature autocorrelation function
In this section, we investigate the effect of fragmentation on infinite-temperature autocorrelation functions under Lindblad evolution, where O ≡ Tr(ρO).The evolution of an operator O is given by O(t) = e tL † (O), which reduces to O(t) = e iHt Oe −iHt without dissipation.
For the observables we consider in the following, the disconnected part is always zero.For closed systems, the infinite-time average of autocorrelation functions is lower bounded by the Mazur bound [72][73][74], which relates a finite saturation value with the presence of conserved quantities.For example, for the family of PF models and considering the local observable O = S z j in a closed system, this bound is given by [46] M PF (S z j ) = where D α is the dimension of Krylov subspace K α .Here {Π α } span a full set of conserved quantities for the Abelian commutant C PF .Ref. [46] numerically found that the bound M PF scales as 1/N in the bulk, hence vanishing in the thermodynamic limit.In Fig. 4, we show the evolution of infinite-temperature autocorrelation functions S z N/2 (t)S z N/2 (0) of the TL model for both closed and open quantum dynamics under different dissipative couplings.For closed systems (green solid line), we numerically evaluate the infinitetemperature correlations by uniformly sampling initial Haar random states as prescribed by quantum typicality [75,76], which saturates to a finite value.For the open dynamics, we simulate the Lindblad evolution using the time-evolving block decimation (TEBD) algorithm [77][78][79], with the infinite temperature configuration as initial state ρ 0 ∝ 1.Under dephasing noise (down-pointing triangles), we find that the autocorrelation function saturates to a lower value than the TL model in closed systems, indicating that the dephasing noise reduces the symmetries of the TL model.The saturation value is exactly the Mazur bound M PF in closed systems given by Eq. ( 16) (blue doted-dashed line).In the inset of Fig. 4, we numerically verify that the saturation values decay as 1/N as previously found in Ref. [46].Appendix A 2 contains additional results for boundary correlations, where a finite saturation value is found.
This agreement between the saturation of autocorrelation of the TL model under dephasing noise and the PF Mazur bound can be explained using the same analysis as for the stationary state ρ ss but now for the stationary value of an operator O(∞) = lim t→∞ e tL † O, which is given by Here O α is a constant given by the overlap of the operator O and the projector.Using O = S z j , we obtain the saturation value of S z j (∞)S z j (0) c as the inner product between S z j and its stationary value S z j (∞) recovering Eq. ( 16).This explains why the autocorrelation function under dephasing noise saturates exactly to the Mazur bound for the PF model.We provide a different proof to the same result in App.B 3, by generalizing the Mazur bound to open systems for diagonalizable L with strong symmetries.

C. Logarithmic negativity
We now investigate the spreading of quantum correlations across the system using the logarithmic negativity [80], an entanglement measure for mixed states defined as Here A 1 = Tr √ A † A is the trace norm, and ρ T B is the partial transpose with respect to a sub-region B, which is given as The logarithmic negativity is an entanglement monotone, which means that it is non-increasing under local quantum operations and classical communication [80], and it is zero for all separable states, i.e., states of the form We study the dynamics of E N starting from the initial state which lies in the largest Krylov subspace associated to . This corresponds to the fully-paired, i.e., trivial dot pattern, subspace.Figure 5 shows the time evolution of the logarithmic negativity.At short times t deph 1/γ, E N increases since the evolution is dominated by the unitary part.However, for t t deph the dephasing noise dominates the dynamics, destroying quantum correlations and leading to a vanishing E N .While our numerical simulations suggest that in the presence of conserved quantities E N has a slow decay, we leave a more detailed analysis for future work.
In fact, the stationary state under dephasing noise, Eq. ( 13), is a separable state for arbitrary initial states.It is the sum of projectors onto product states in the local z basis |ψ αβ = |ψ A αβ ⊗ |ψ B αβ , appearing as a result of the classical fragmentation and Hermitian jump operators.Hence, it can be written as αβ |.Therefore, the logarithmic negativity for an arbitrary bipartition with an arbitrary initial state is zero.This result generalizes to stationary states for systems with Abelian commutants spanned by a local product basis and Hermitian jump operators.

D. Operator space entanglement
While quantum correlations eventually vanish in the presence of dephasing noise, information continues its spreading in the presence of conserved quantities.We characterize this spreading using the operator space entanglement (OSE), which measures the von Neumann entropy of the vectorized density operator ρ → |ψ(ρ) , using Choi's isomorphism |σ i σ i | → |σ i σ i [81].With the Schmidt decomposition of |ψ(ρ) , the OSE is given by, where the Schmidt values λ a are normalized to a λ 2 a = 1.In the presence of conserved quantities, the OSE can be split into two types of entanglement: the number entanglement S num and the symmetry-resolved entanglement S res [82,83], S num is the Shannon entropy associated with the fluctuations of the conserved quantities in half of the system and S res the weighted von Neumann entanglement entropy within each symmetry sector.
We study the evolution of OSE starting from the same initial state |ψ 0 = ⊗ j |+ j in Fig. 6.For small γ = 0.1, similarly to the logarithmic negativity, the OSE grows for a time t 1/γ, and is then suppressed by the dissipation.However, the OSE saturates to size-dependent finite values (Fig. 6a).For large γ = 10, the OSE is largely suppressed, which allows for efficient TEBD simulation for larger system sizes.We observe that the OSE grows even with the presence of dissipation until saturation (Fig. 6b).The saturation values can be calculated from the expression for the stationary state ρ ss in Eq. (13).Vectorizing the stationary density matrix ρ ss → |ψ ss , one finds that the saturation value of the OSE is given by the von Neumann entropy of the state |ψ ss , which was analytically obtained in Ref. [69].In particular, it was shown that S OP (ρ ss ) = S num (|ψ ss ) scales as O( √ N ) with system size N and that S res (|ψ ss ) = 0 (Fig. 6c).
A recent study argued that the OSE grows logarithmically in the presence of a U(1) charge validating this expectation for certain systems [83].For the U(1)conserving XXZ chain considered in Ref. [83], the authors found that the strongly dephased dynamics can be approximated by a symmetric simple exclusion process of hardcore particles.There particle fluctuations across the bipartition resulted in a logarithmic growth of the number entropy, while the symmetry-resolved entanglement vanished.
In the following, we extend this analysis to the presence of the non-local conserved quantities that characterize the fragmented structure of the stationary state, which helps to understand the OSE growth observed in Fig. 6b.Unlike Ref. [83], the number entropy of the systems we consider in this work is related to the fluctuations of non-local conserved quantities, the color-dot patterns.Analogously to the U(1) charge N c that admit the decomposition N c = N L + N R , we split the global dot pattern A k into left and right patterns such that we can keep track of their fluctuations.For example, for the fully-paired state | , the left and right dot patterns after a half-chain bipartition are given by A k = { } and Āk = { }, respectively.Similar to the case of zero total charge with N R = −N L , the right dot pattern is a reflection of the left dot pattern for the fully-paired subspace.As a result, the half-chain number entanglement entropy is given by with p A k the probability of having the left dot pattern A k3 .Operator space entanglement and number entanglement under dephasing noise.The initial state is |ψ0 = ⊗j|+ j .(a) Lindblad dynamics of the OSE under dephasing noise Lj = S z j with γ = 0.1 using ED.For small γ, the OSE increases at short times t 1/γ when the dynamics is governed by the unitary term, then decreases and saturates to a size-dependent value.(b) Lindblad dynamics with large γ = 10 using TEBD.The OSE is largely suppressed by the dissipation, which allows efficient TEBD simulation.The data suggests a logarithmic growth with a rate increasing over time (see main text).(c) The analytic results of the OSE for the stationary state (black dots), the saturation values of SOP under Lindblad dynamics (upper-pointing triangles), and Snum under stochastic dynamics (down-pointing triangles) show quantitative agreement.The saturation values under Lindblad dynamics are obtained with the same TEBD parameters as in (b).The OSE of the stationary state in Eq. ( 13) scales as O( √ N ) with system size.(d) Number entanglement of the effective stochastic dynamics, which shows similar behavior as in the Lindblad dynamics with large γ.Each curve is averaged over 10000 random samples.
In the limit of strong dephasing, we can derive an effective Lindblad evolution using degenerate perturbation theory for open quantum systems [84].We do so by splitting L = L 0 + L 1 into the unperturbed contribution L 0 and the perturbation L 1 in the limit |J j |/γ → 0. Here, Since the initial state |ψ 0 = ⊗ j |+ j lies in the fully-paired subspace of the PF model, the stationary states of L 0 are given by ρ σ 0 = |σ σ|, where |σ are all possible fully-paired product states.The perturbation L 1 breaks this degeneracy inducing transitions among different ρ σ 0 .Performing the perturbation theory to second order in |J j |/γ we find the effective Liouvillian [47,[83][84][85]] where P is the projection onto the subspace spanned by ρ σ 0 .This effective dynamics reduces to a classical Markov evolution ∂ t ρ(t) = −W eff ρ(t) for the diagonal components of ρ in the fully-paired product basis ρ σ 0 with W eff = j g αβ j (|αα ββ|) j,j+1 is the Markov generator given by a PF model with coefficients g αβ j obtained in App.C 1.This implies that the effective dynamics indeed preserves the commutant algebra associated to the PF model C PF .
For an XXZ model under dephasing noise in Ref. [83], the corresponding effective stochastic evolution can be mapped to a simple exclusion process, from where an analytical prediction for the growth of S num could be obtained.However, we are not aware of any analysis of the evolution generated by W eff .Hence, we numerically simulate it in a manner that can be compared to the implementation for open quantum dynamics.In the basis spanned by {ρ σ 0 }, the probability vector with entries p σ (t) at discrete time t is given by p σ (t) = σ (P t ) σσ p σ (0) where P = e −W eff [47].Transition probabilities are given by the corresponding entry in the matrix P , which is symmetric, and satisfies P σσ ∈ [0, 1] together with σ P σσ = 1.Hence, detailed balance holds with respect to a stationary state, which is the uniform distribution over all fully-paired states.This corresponds to the stationary state ρ ss of the Lindblad dynamics.To efficiently implement this evolution, we consider a brick-wall circuit structure where 2-site local gates P j,j+1 randomly permute among two-site local spin configurations in the z-basis as, e.g., in Refs.[27][28][29][30][31][32]44].Starting from the initial product state ⊗ j |+ j , we then compute the evolution of the number entropy S num as given in Eq. ( 22) by averaging over various circuit realizations.More details about the numerical implementation can be found in App.C 2.
In Fig. 6, we compare the open quantum dynamics (panel b with γ = 10) with the stochastic one in panel d.The latter allows us to simulate larger system sizes and longer times than what is accessible by TEBD simulations.We observe a growth of the number entanglement of the stochastic model in Fig. 6b, which agrees with the numerical results obtained in the quantum setup.However, we are unable to provide an analytical prediction for the observed scaling of growth as for the U(1)-symmetric systems.Assuming a logarithmic growth of the OSE S(t) = S 0 + η log(t), we find that the growth rate η slightly increases over time.Note that a similar effect is also observed in Fig. 2a of Ref. [83] for U(1) symmetric systems, which is caused by finite time effects.Our numerical simulations reach a saturation value for the S num (red down-pointing triangles) that agrees with the analytical result (black dots) and the saturation of the OSE under the quantum Lindblad dynamics (blue upper-pointing triangles) as shown in Fig. 6c.

V. STRUCTURE-PRESERVING NOISE
In the previous section, we observed that the dephasing noise reduced the QF of the TL model to the classical one.This led to vanishing quantum correlations as measured by the E N , while classical correlations (S num ) could still propagate due to fluctuations of the remaining conserved quantities.We now consider a dissipative bath preserving the QF and investigate the effects of the system being fragmented in an entangled basis.We choose L j = e j,j+1 acting on two consecutive sites, which is an element of the bond algebra A TL .Hence, the open commutant algebra agrees with that of the TL model

A. Stationary states with quantum fragmentation
When considering quantum structure-preserving noise, the stationary state inherits the QF of the TL model leading to the general expression where (M λ ) αα = Tr(Π λ α α ρ 0 ) is the d λ × d λ matrix of overlaps between the initial state ρ 0 and Π λ αα with Π λ α ≡ Π λ αα .There are two major differences which distinguishes this from the stationary state discussed in the previous section.First, there are stationary phase coherences, i.e., L(Π λ αα 1) = 0, captured by the non-zero overlaps with the conserved intertwine operators Π λ αα .Recall that these appear as a consequence of C TL being non-Abelian.As in the case of dephasing noise, the conserved projectors give the stationary state Π λ α 1/D λ in the diagonal subspaces.These projected identities indicate full decoherence within the subspaces H (A) λ induced by L j ∈ A. Nonetheless, intertwine operators acting on the off-diagonal subspaces, guarantee non-vanishing coherences for generic initial states [67], indicating that the whole system does not fully decohere.Figure 7 shows an example of the Lindblad evolution for the initial state in Eq. ( 14) displaying non-zero overlap onto the nondegenerate fully-dimerized subspace (λ = 0) and onto two degenerate Krylov subspaces (λ = 1) in the entangled basis of the TL model.The system evolves to the stationary state with projected identities both in the diagonal and off-diagonal degenerate subspaces.Second, the projected identity Π αα 1 within each Krylov subspace is a mixture of entangled basis states.This implies that the stationary state is typically not separable unless for fine-tuned initial states.As we find in the following, this is also signalled by the behavior of the logarithmic negativity.Moreover, the exponentially large (in system size) dimension of the commutant algebra as a consequence of HSF, dim(C) = λ d 2 λ ∼ e aN , turns into a strongcoherent in the case of non-Abelian Cmemory of the initial configuration.Information about the initial state is stored by the weight on the invariant subspaces λ are decoherence-free subspaces and noiseless subsystems immune to dissipation, which are extensively studied in the context of error correction and fault tolerant quantum computation [86][87][88][89][90][91].

B. Infinite-temperature autocorrelation function
Once again we can use a similar analysis to that of the stationary state to derive the saturation value of the spinspin autocorrelation function S z j (t)S z j (0) .One finds that a general operator O relaxes to the stationary value where (O λ ) αα = Tr(Π λ α α O).Here, O λ is a d λ × d λ matrix with elements given by the overlap of the operator and the corresponding projector or intertwine operator.
Therefore, for a local operator S z j , the saturation value of the autocorrelation is given by Tr(S z j (∞)S z j (0))/3 N , which is exactly the Mazur bound of the TL model for unitary evolution This agrees with the numerical results shown in Fig. 4, where the autocorrelation functions saturate to the same value for the closed system (green solid line) and under the structure-preserving noise (upper-pointing triangles).The finite-size scaling of the saturation values suggests that it is not vanishing either in the bulk or at the edge (see App.A 2).

C. Logarithmic negativity and operator space entanglement
A vanishing or non-vanishing bulk autocorrelation function is not sufficient to distinguish classical from quantum fragmentation.For example, the bulk autocorrelation functions decay to zero for the t − J z chain but remain finite for certain dipole-conserving models, both of which are classical fragmented [16,19].
However, a sharp contrast can be detected in the behavior of the logarithmic negativity in the presence of different types of baths.While we found a vanishing negativity for dephasing noise when starting from the initial state ⊗ j |+ j , we find that E N saturates to a size-dependent value at long times under the structurepreserving noise, indicating that the system evolves towards an entangled stationary state (see Fig. 8).Moreover, the scaling of the negativity with system size is directly computed from the stationary state in Eq. ( 26) and shown in the inset of Fig. 8, suggesting that the stationary state satisfies a volume-law.The source of this non-vanishing value is the fact that the system is fragmented in an entangled basis, hence providing a clear signature to distinguish quantum and classical fragmentation.
Thus we propose the logarithmic negativity of stationary states as a probe to distinguish quantum from classical fragmentation.Generally, identifying CF structure is an easier task that can be achieved by iteratively applying local terms of Hamiltonian to a root product state.However, there can still be a finer structure within these Krylov subspaces due to quantum fragmentation appearing in an unknown entangled basis.To detect whether such a finer structure exists, one could start from an initial state within a Krylov subspace, and study the dynamics of the logarithmic negativity under a dissipative bath, which should preserve all the symmetries of the Hamiltonian.This means that the jump operators should be elements of the bond algebra L j ∈ A and Hermitian.While systems showcasing only CF evolve towards a separable stationary state with zero negativity, systems that are quantum fragmented can lead to non-zero logarithmic negativity.
Before concluding this section, we study the evolution of the OSE and compare its saturation value to that obtained from the stationary state.The results are shown in Fig. 9 for γ = 0.1 (panel a) and γ = 10 (panel b).In this case, the dynamics of OSE cannot be efficiently studied even in the regime γ 1 for the following reasons.Logarithmic negativity under structurepreserving noise.Time evolution of the logarithmic negativity using ED under the structure-preserving noise Lj = ej,j+1, with γ = 0.1.The initial state is |ψ0 = ⊗j|+ j .The logarithmic negativity EN increases at t 1/γ and then saturates to a finite size-dependent value.In the inset, we show the scaling of this saturation value as directly computed from the stationary state in Eq. ( 26), also included in the main panel (dashed green lines).While our numerical simulations are limited to system sizes N ≤ 10, the scaling with system size suggests volume-law.
First, the stationary-state subspace of the unperturbed contribution L 0 is spanned by entangled states.Obtaining an orthonormal set of these entangled states requires full diagonalization of L 0 .Second, transition among entangled states cannot be modelled by local updates on local configurations and hence cannot be mapped to a classical stochastic circuit evolution.Moreover, while under dephasing noise we could directly extract the number entanglement by calculating the probabilities of the dot patterns of the product states, this is not the case for entangled states which involved entangled dot patterns.This raises the general question whether one can capture quantum fragmentation phenomena using classical stochastic dynamics.

VI. CONCLUSIONS AND OUTLOOK
The goal of our work was to examine how HSF impacts open Lindblad dynamics, taking into account whether the coupling to the bath maintains or disrupts fragmentation in an entangled basis.By analyzing the symmetries of the Liouvillian, we were able to analytically derive the stationary state and characterize the dynamics of autocorrelation functions and entanglement combining analytical and numerical methods.First, we found that for a dephasing noise -that reduces the quantum fragmentation of the TL model to the classical fragmentation of the PF-the stationary state is a separable state with zero quantum correlations.This holds generically for classically fragmented open systems with Hermitian jump operators.Nonetheless, the OSE increases as a function of time due to the fluctuations of the non-local conserved charges as captured by an effective stochastic evolution in the regime of strong dephasing.On the other hand, for a dissipative coupling preserving the QF of the TL model, the system evolves to a highly-entangled stationary state with size-dependent logarithmic negativity.This finite saturation value is a dynamical property distinguishing classical from quantum fragmentation in open quantum systems, while for unitary evolution both classical and quantum fragmentation lead to volume-law entanglement entropies.In addition, there exist stationary coherences in the off-diagonal subspaces due to non-Abelian commutant algebras, indicating that the system does not fully decohere.Although the system shows distinct entanglement properties under the two couplings, finite autocorrelation functions could persist under both types of dissipation.Moreover, the extensive fragmentation of the Hilbert space translates into exponentially many (in the volume of the system) degenerate stationary states signaling a strong dependence on the initial state.
The preceding discussion has highlighted three critical components: (1) the distinction between classical and quantum fragmentation, which is synonymous with product and entangled basis spanning the fragmented structure respectively.This translates into stationary identity matrices within Krylov subspaces in terms of either product or entangled states, respectively, where the latter leads to a finite negativity at long times.(2) The distinction between Abelian and non-Abelian commutants; A non-Abelian commutant results in the presence of stationary coherences, which indicates a coherent memory of the initial state [52,92]; and (3) the exponential dimension of the commutant as caused by HSF, which leads to a large degeneracy of stationary states and a strong dependence on the initial state.
For future work, it will be interesting to understand whether similar entanglement dynamics as the one found for quantum fragmented systems, appears for polynomially large commutants.For example, conventional symmetries such as SU(2) also lead to a decomposition of the Hilbert space into symmetry sectors spanned by an entangled basis, which may evolve to a stationary state with finite negativity for specific initial states.However, with exponentially large subspaces that scale as the size of the Hilbert space, the stationary state is highly mixed, which can exhibit a different dependence of entanglement with system size.
We also leave it open to explore classical and quantum fragmentation in the presence of weak symmetries [52,53].In fact, an example of classical (local) fragmentation in this weak sense already appeared in Ref. [51].A natural adaptation of the commutant algebra formalism consists of considering the vectorized form of the Lindbladian L → L acting on the Hilbert space H⊗H and define the commutant as the set of (super)-operators commuting with every local term of L. For example, it would be interesting to understand whether there are examples of quantum fragmentation and non-Abelian commutants for weak symmetries, and if so, whether they lead to similar phenomenology as the one found in this work.
Finally, while several recent studies [27][28][29][30][31][32][43][44][45]] have employed block (local) cellular automaton dynamics to investigate the impact of classical fragmentation on infinite-temperature correlations, our work raises the following question: Is it possible to construct a blocked cellular automaton with finite-size gates that simulates the dynamics and capture the entanglement properties of quantum fragmentation?If it is not possible, the obstruction to find such cellular automaton could be used as a definition of quantum fragmentation.

ACKNOWLEDGMENTS
The tensor-network calculations in this work were performed using the TeNPy Library [93] Data and materials availability.Data analysis and simulation codes are available on Zenodo upon reasonable request [94].The TL model exhibits QF (in an entangled basis) where the Krylov subspaces are labeled by product or entangled dot patterns.In addition, due to the non-Abelian commutant algebra C TL , the Krylov subspaces with dot patterns of the same length are degenerate.
We provide some simple examples of how to construct the entangled basis of the TL model.We label the basis states by |ψ λ αβ , where 2λ is the number of dots, α = 1, ..., d λ denotes different degenerate Krylov subspaces for fixed λ, and β denotes different basis states in the same Krylov subspace.For a system with two sites N = 2, the fully-dimerized Krylov subspace with λ = 0 (zero dots) is one dimensional, with |ψ 0 1,1 = | .For λ = 1 with two dots, the Krylov subspaces are also one-dimensional with |ψ 1 α,1 = | , such that e j,j+1 | = 0.The dot state can be a product state, |σ 1 σ 2 with σ 1 = σ 2 , or an entangled state such as 1 √ 2 (| + + − | − − ).The Krylov subspaces with N = 2 and λ = 1 have a degeneracy of d 1 = 8, i.e., there are in total eight different dot patterns which consist of two dots.Note that the choice of dot patterns is not unique, any linear superposition of dot patterns works.For larger system sizes with Krylov subspaces of dimension D λ ≥ 1, we apply e j,j+1 on a root state of the subspace to generate other basis states.For example, for N = 4 with λ = 1, the Krylov subspace is threedimensional with basis states The dot pattern ( ) is conserved and labels this Krylov subspace.A systematic way to construct the complete basis is given by Ref. [71].

Finite-size scaling of the autocorrelation functions
Both the PF and TL models exhibit strong fragmentation [16,46].Figure 10a shows that the number of Krylov subspaces for the PF and the TL models scales exponentially with the system size [69,71].Figure 10b shows that the ratio between the dimension of their largest Krylov subspace and the total Hilbert space dimension scales as D max /D ∼ exp(−aN ).We study the non-ergodic behavior due to strong fragmentation with the long-time average of autocorrelation functions, which is given by We study C z j (∞) with random circuits using ED, which is shown in Fig. 10(c-d).At the boundary, C z 0 (∞) decays with the system size for both TL and PF model, but saturates to a finite value in the thermodynamic limit.This indicates that there are infinite coherence times at the boundary for both classical and quantum fragmentation.At the bulk, as discussed in the main text and in Fig. 4, the autocorrelation functions C z N/2 (∞) of the PF and TL models coincide with the saturation values in open systems under dephasing noise and the structure-preserving noise, respectively.The bulk autocorrelation decays as 1/N and vanishes for the PF model, while for the TL model, the numerical results suggest that the autocorrelation functions saturate to finite values.

Derivation of non-equilibrium stationary states
With the analysis of the commutant algebra of the Lindblad system, we obtained a full set of conserved projectors Π α , which decompose the operator space into minimal subspaces with independent dynamics.Now we prove the uniqueness of the eigenstates with zero eigenvalues (fixed points) of the Liouvillian within the minimal subspaces.This can be explained as follows [67,68]: Density matrices form a convex set S, where the boundary ∂S consists of all states with a lower rank [95].Assume that both ρ 1 ∞ and ρ 2 ∞ are stationary states in one diagonal minimal block B α .Due to the linearity of the Lindblad equation, ∞ forms a line of stationary states.Assume that the line intersects with the boundary ∂S at ρ 3 ∞ , which has rank(ρ 3 ∞ ) smaller than the dimension of the subspace.The range of ρ 3 ∞ is then a smaller subspace that contains a stationary state.This indicates that we can further decompose B α , which is a contradiction with the fact that B α is a minimal subspace.Therefore, within each diagonal minimal block, there is at most one stationary state.In our case, we have proven that the unique stationary state within the subspace is the projected identity, Π α 1/D α ∈ B α .
In the off-diagonal subspaces, the existence of fixed points, i.e., the stationary coherences, is guaranteed by the conserved intertwined operators Π λ αα between two degenerate diagonal subspaces for non-Abelian commutant algebras.This is given in Theorem 18 of Ref. [67].There is also a unique fixed point in each off-diagonal subspace.Otherwise, the conserved intertwine operator gives an extra fixed point in the diagonal subspace, which is a contradiction.
To obtain the general expression of the stationary state in the full operator space, we perform a spectral decom-position of the Liouvillian superoperator.Due to the non-Hermiticity of L, there is in general a different set of eigenstates for L and L † given by [52,65] or equivalently (σ † n )L = λ n σ † n for the latter, i.e., different left and right eigenstates of L. They satisfy the biorthonormal relation σ m |ρ n ≡ Tr(σ † m ρ n ) = δ mn .The left and right eigenmatrices span a full basis, such that we can expand the initial state as ρ 0 = n c n ρ n , with c n = Tr(σ † n ρ 0 ).The eigenspectrum of L consists of eigenvalues with Re(λ) ≤ 0. Consider systems without purely-imaginary eigenvalues: in the long time limit, the dynamics is governed by eigenmatrices with zero eigenvalues.Therefore, the full stationary state is then given by [52] With the analysis of the Hilbert space structure, we have identified the full set of eigenstates with zero eigenvalues of L, which are the stationary states and stationary coherences {Π αα 1}, as well as their corresponding conserved quantities {Π αα }.With Eq. (B2), we obtain the general expression for the stationary state of fragmented systems, specified in the main text by Eq. ( 13) and Eq. ( 26) in the case of classical and quantum fragmentation respectively.

Mazur bound in open system
The Mazur bound in closed systems relates the infinitetime average value of autocorrelation functions (nonergodicity) to the presence of conserved quantities [72][73][74] with D as the dimension of Hilbert space, and Y = Y † .With strong symmetries, the set of conserved quantities {J µ } satisfy [H, J µ ] = [L j , J µ ] = 0 for all L j , indicating that L|J µ = 0. Therefore, the corresponding right eigenmatrices are also given by J µ4 .In the long-time average, all the oscillating (Imλ = 0) and decaying terms (Reλ < 0) vanish.Therefore, the contributions will be given by the {J µ } which are associated with zero eigenvalues, where we define Y (t) ≡ lim T →∞ 1 T T 0 Y (t)dt.This is the Mazur bound in the open system, with a set of orthogonal conserved quantities {J µ }.In general, for conserved quantities {Q µ } not orthogonal, the Mazur bound can be written as where (K) µν ≡ Q † µ Q ν is the correlation matrix.For open fragmented systems, with the choice of L j ∈ A, all elements in C commute with the Hamiltonian and jump operators.A full set of orthogonal conserved quantities in C are the projectors onto the Krylov subspaces {Π α } and the intertwine operators {Π αα }.With jump operators S z j and e j,j+1 , the Mazur bound gives the same results derived from the stationary states in Eq. ( 16) and Eq. ( 28), respectively.

Simulation of the stochastic dynamics
The effective dynamics of the TL model under dephasing noise can be mapped to a classical Markov process, with the stochastic matrix W eff in Eq. (C2).Here we provide more details of our numerical simulation.In correspondence to the random circuit setting for the quantum Lindblad dynamics, we implement the dynamics by classical circuits, with two-site gates U t,j permuting among the classical configurations σ = {s 1 , ..., s N }, where s i ∈ {+, 0, −}.These classical configurations σ correspond to ρ σ 0 with |σ = |σ 1 , ..., σ N .A two-site gate U t,j acting on the configuration σ = {..., s j , s j+1 , ...} at time t gives a new configuration σ = {..., s j , s j+1 , ...} with a transition probability.A configuration with s j = s j+1 can transform to a new configuration with s j = s j+1 , with the transition probability given by the probability matrix s j s j+1 |P j,j+1 |s j s j+1 with P j,j+1 = exp(− J 2 j γ M i,i+1 ) [47].This is a pair-flip action.For s j = s j+1 , the configuration is unchanged.To compare with the Lindblad dynamics, we start from the initial configuration with all s i = +.Averaging over random circuit realizations, we obtain the time evolution of the probability of dot patterns p A k (t), and thus the number entanglement S num (t).The mapping of the effective dynamics to the stochastic dynamics allows simulation for much larger system sizes and longer times.The half-chain entanglement of |ψ ss was studied in Ref. [69].It was shown that the symmetry-resolved entanglement S res = A k p A k S res (A k ) = 0, as after resolving the left dot pattern A k , all the configurations contribute equally (thus the state can be written as a product state with S res (A k ) = 0 for all A k .Hence, only the number entropy given by S = − A k p A k log p A k , with A k as the dot pattern of the left part of the chain remains.For large system sizes N , the entanglement scales as S ∼ O( √ N ).In Fig. 11, we show the Lindblad evolution of S num and S res under dephasing noise for γ = 10.We show that the symmetry-resolved entanglement is small compared to the number entanglement during time evolution, and saturates to zero for the stationary states.
, we denote the spins with different colors as | = |+ , | = |0 , and | = |− .Using this representation, the pair-flip terms of H PF change neighboring spins with the same color,

B 2 .FIG. 3 .
FIG.3.Time evolution of the density matrix under dephasing noise.Time evolution of the density operator under dephasing noise Lj = S z j using ED, with system size N = 4 and γ = 1.The initial state is specified in Eq. (14), having nonzero overlap with three different Krylov subspaces.The color intensity in the figures is the magnitude of matrix elements |ρij| in the product z basis, and the solid blue lines separate different Krylov subspaces.
FIG. 4.Infinite-temperature autocorrelation functions.Time evolution of bulk autocorrelation functions S z N/2 (t)S z N/2 (0) for unitary and open quantum dynamics using different jump operators Lj.We use a system size N = 12 and γ = 1.The unitary dynamics is calculated using ED, while we use TEBD for the Lindblad dynamics with bond dimension χ = 256.Under the noise Lj = ej,j+1 that preserves the quantum fragmentation, the autocorrelation function saturates to the same value as in the closed system for the TL model HTL.Under the dephasing noise Lj = S z j , the autocorrelation function saturates to a finite value, which is the Mazur bound of the PF model (blue dotted-dashed line).The spin-flip noise Lj = S x j further destroys the classical fragmentation, which leads to vanishing autocorrelation functions.

8 FIG. 5 .
FIG.5.Logarithmic negativity EN under dephasing noise.Time evolution of the logarithmic negativity as given in Eq. (18) using ED, under dephasing noise Lj = S z j with γ = 0.1.EN increases at short times t 1/γ, when the dephasing noise kicks in destroying quantum correlations.
FIG. 6.Operator space entanglement and number entanglement under dephasing noise.The initial state is |ψ0 = ⊗j|+ j .(a) Lindblad dynamics of the OSE under dephasing noise Lj = S z j with γ = 0.1 using ED.For small γ, the OSE increases at short times t 1/γ when the dynamics is governed by the unitary term, then decreases and saturates to a size-dependent value.(b) Lindblad dynamics with large γ = 10 using TEBD.The OSE is largely suppressed by the dissipation, which allows efficient TEBD simulation.The data suggests a logarithmic growth with a rate increasing over time (see main text).(c) The analytic results of the OSE for the stationary state (black dots), the saturation values of SOP under Lindblad dynamics (upper-pointing triangles), and Snum under stochastic dynamics (down-pointing triangles) show quantitative agreement.The saturation values under Lindblad dynamics are obtained with the same TEBD parameters as in (b).The OSE of the stationary state in Eq. (13) scales as O( √ N ) with system size.(d) Number entanglement of the effective stochastic dynamics, which shows similar behavior as in the Lindblad dynamics with large γ.Each curve is averaged over 10000 random samples.

>FIG. 7 .
FIG. 7. Time evolution of the density matrix under structure-preserving noise.Time evolution of the density operator under structure-preserving noise Lj = ej,j+1 using ED, with system size N = 4 and γ = 1.The color indicates the matrix elements |ρij| in the entangled basis.The stationary state consists of projected identities in all diagonal blocks and degenerate off-diagonal blocks, while all elements in non-degenerate blocks vanish.
FIG. 8.Logarithmic negativity under structurepreserving noise.Time evolution of the logarithmic negativity using ED under the structure-preserving noise Lj = ej,j+1, with γ = 0.1.The initial state is |ψ0 = ⊗j|+ j .The logarithmic negativity EN increases at t 1/γ and then saturates to a finite size-dependent value.In the inset, we show the scaling of this saturation value as directly computed from the stationary state in Eq. (26), also included in the main panel (dashed green lines).While our numerical simulations are limited to system sizes N ≤ 10, the scaling with system size suggests volume-law.

FIG. 9 .
FIG. 9. OSE under structure-preserving noise.Time evolution of the OSE under the structure-preserving noise Lj = ej,j+1 with the initial state |ψ0 = ⊗j|+ j .Data is obtained using ED for γ = 0.1 and TEBD for γ = 10 with maximal bond dimension χ = 1000.(a) For γ = 0.1, the OSE increases at short times and then decreases and saturates to a finite value.The inset shows the OSE of the stationary state obtained from Eq. (26) (circle), which matches the saturation values obtained by ED (upper-pointing triangle).(b) For γ = 10, the OSE saturates to the same values as for γ = 0.1.
Appendix A: Fragmentation of PF model and TL model 1.Entangled fragmentation basis of the TL model

FIG. 10 .
FIG. 10.Strong fragmentation and saturation of autocorrelation functions.(a) Number of Krylov subspaces K = λ d λ for PF and TL model, which scales exponentially with the system size.(b) Ratio between the dimension of the largest Krylov subspace Dmax and the dimension of the total Hilbert space D = 3 N , which scales as Dmax/D ∼ exp(−aN ) with a < 1.This indicates that both models exhibit strong fragmentation.(c-d)Finite-size scaling of the long-time average of the autocorrelation functions at the boundary (j = 0) and at the bulk (j = N/2) as 1/N .At the boundary, both the PF and TL model show infinite coherence time in the thermodynamic limit.At the bulk, the TL model has nonvanishing autocorrelation functions, while for the PF model, the autocorrelation function decays as 1/N and vanishes.

FIG. 11 .
FIG. 11.Number entanglement and symmetryresolved entanglement.Time evolution of the number entanglement Snum (square) and the symmetry-resolved entanglement Sres (circle) for different system sizes N by ED.The Sres is small compared to Snum, and saturates to zero.
. The authors thank Barbara Kraus, Olexei Motrunich and Sanjay Moudgalya for valuable advice.Y. L. thanks Fabian Essler, Zongping Gong, Johannes Hauschild, Dieter Jaksch, Yujie Liu, Hannes Pichler, Elisabeth Wybo, and Zhongda Zeng for helpful discussions.P.S. acknowledges support by the Walter Burke Institute for Theoretical Physics at Caltech, and the Institute for Quantum Information and Matter.This research was financially supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 771537.F.P. acknowledges the support of the Deutsche Forschungsgemein-schaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC-2111-390814868.F.P.'s research is part of the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus.
Here we prove that the dynamics of the TL model under dephasing noise characterized byC O deph = {e i,i+1}, {S z j } is exactly C PF .First, e i,i+1 and S z j are elements of A PF , therefore, A O deph ⊆ A PF .Second, all elements of the PF algebra can be generated by the local terms of the TL model and the dephasing noise.].The linear combinations of such products of S z j and e i,i+1 produce all the local terms of the PF model, which indicates that A PF ⊆ A O deph .Hence, we have A O deph = A PF , which also means that they have the same commutant, C O deph = C PF .Altogether one finds that the fragmentation structure of the TL model under dephasing noise is determined by C PF .
. In the main text we derived the saturation value of autocorrelation function O(∞)O(0) by the stationary value of the operator O, which coincides with the Mazur bound of the closed systems.The same conclusion can be achieved by generalizing the Mazur bound to open systems.
i } and {r i } form a complete set of basis, i.e. i |r i l i | = 0 and satisfy the biorthonormal relation.Thus we can expand arbitrary observable Y as |Y = i |r i l i |Y [65].The autocorrelation function is