On the generality of symmetry breaking and dissipative freezing in quantum trajectories

,


Introduction
Symmetries are fundamental to our understanding of nature.Quantum physics is no exception to this, and the invariance of the Hamiltonian of a closed quantum system under a given transformation can have far-reaching implications, which include eigenstate degeneracy [1], an absence of thermalisation [2,3] and the protection of topological order [4,5].
In an open quantum system, the influence of an external environment complicates matters and must be accounted for when identifying symmetries and assessing their impact on the system.The last decade has seen significant progress in this regard, with extensive classifications of the symmetries of open systems having been achieved [6][7][8][9][10] alongside an understanding of their physical consequences in terms of steady state degeneracy [11,12], transport properties [13][14][15], non-stationarity [10,16], quantum synchronisation [17,18], off-diagonal long-range order [19,20] and dissipative phase transitions [21][22][23].Results such as this have made it clear that the combination of symmetries and an external environment provides a new pathway for realising unique, exploitable states of matter in a quantum system.
Recently, a new symmetry-based phenomenon, termed 'dissipative freezing', was witnessed in a large quantum spin precessing under the influence of dissipation parallel to the axis of rotation [24].It was observed how the individual quantum trajectories -solutions to the stochastic Monte Carlo unravelling of the master equation -randomly selected a specific symmetry sector to converge into in the long-time limit; breaking that symmetry at the trajectory level.Due to the simplicity of the model at hand, an analytical proof of this unexpected result was provided, although it was unclear whether this phenomenon could be observed in more general setups.Since this result, several independent works have noted the manifestation of this freezing in other open system setups [25,26], but without offering a deeper explanation for its origin.
In this work, we consider the question of whether dissipative freezing emerges as a general feature of open quantum systems with a strong symmetry.We develop several, simple, mathematical and physical perspectives on the problem to reasonably justify that this is the case, as well as uncovering the direct connection between this phenomenon and previous, rigorous, mathematical proofs establishing the ergodicity of open systems when the external environment is subject to continuous monitoring.We follow these insights up with the derivation of a lower-bound for the rate-of-freezing in terms of the spectral properties of the Liouvillian.Finally, we give several illustrative numerical examples from a range of open systems, where dissipative freezing and the insights developed in this paper can be observed.
Importantly, the simple understanding of dissipative freezing provided in this work allows us to identify the exceptional situations in which, despite the presence of a strong symmetry, dissipative freezing is absent.The most prominent of these situations involves the closing of the spectral gap in the off-diagonal symmetry subspaces and heralds the emergence of inter-sector, traceless long-lived modes in the Liouvillian's spectrum.These modes ensure the preservation of coherences and information between symmetry sectors despite the environmental influence and can lead to complex physical phenomena such as non-stationarity and synchronisation [10,17,18].By framing this gap closing in terms of an absence of dissipative freezing, we provide a novel, computationally efficient method for identifying the presence of such intersector traceless modes in open systems.

The GSKL Equation and the Strong Symmetry
Our starting point is a quantum system coupled to an environment under the Markov approximation.The most general completely positive trace-preserving (CPTP) map for the evolution of this system's density matrix ρ(t) at time t is given by the Gorini-Sudarshan-Kosakowski-Lindblad (GSKL) equation [27,28]: where H is the system's Hamiltonian and L 1 , L 2 , ..., L M are a series of 'jump' operators which describe the interaction of the system with the environment.We set ħ h = 1 throughout and use γ j to set the dissipation strength.
Due to the addition of the dissipative term in the equation of motion of the system, symmetries can no longer simply be identified by operators which commute with the Hamiltonian H and there are, in fact, several different types of symmetry that an open system can possess [6].
Here, we focus on the 'strong symmetry', which can be identified by a Hermitian operator A that satisfies ensuring the observable 〈A〉 is a constant of motion.
As A is Hermitian, we can diagonalise it where the index α runs over the D distinct eigenvalues λ α of A and β runs over the corresponding d α degenerate eigenvectors |α, β〉 for a given α, i.e.A|α, β〉 = λ α |α, β〉.Throughout this work we will use α to notate a 'subspace' or 'symmetry subspace' spanned by these degenerate vectors, and define a 'block' as the operator space spanned by the We refer to a given block as 'diagonal' if α = α ′ and 'offdiagonal' otherwise.Now, as the Hamiltonian and the jump operators all commute with A it immediately follows that A(O |α, β〉) = λ α (O |α, β〉) for O = H, L j , L † j ∀ j and so they cannot map an eigenstate of A out of its subspace.They are thus all 'block-diagonal', meaning we can write them in the following form The presence of a strong symmetry therefore immediately implies there exists a simultaneous block-decomposition of the Hamiltonian and jump operators.We emphasize that the converse is also true: if there exists some basis {|α, β〉} in which the Hamiltonian and jump operators are simultaneously block-diagonal then there exists a strong symmetry of the system given by the operator A in Eq. (3) -i.e. the operator which is just some scalar multiple of the identity matrix in each diagonal block and zero everywhere else.Due to this strong symmetry there will be multiple steady states of the system, i.e. multiple density matrices ρ ∞ which satisfy Lρ ∞ = 0 [6].In fact, the basis of steady states (the nullspace of the Liouvillian) will contain at least D trace unity Hermitian steady states.This is because the projection of the Liouvillian into the superoperator space formed from two copies of a given diagonal block constitutes a valid CPTP map and that via Evans' theorem there will exist at least one trace unity steady state for such a map [29,30].There can also exist traceless steady states, or even traceless imaginary eigenmodes (modes ρ imag with imaginary eigenvalues -i.e.Lρ imag = iλρ imag , λ ∈ ̸ =0 ) of the Liouvillian.In the case that these traceless modes are the eigenvectors of the projection of the Liouvillian into the superoperator subspace formed from two copies of a given off-diagonal block, we will refer them as 'intersector traceless, non-decaying modes'.This is because they are solely comprised of coherences between states in different symmetry sectors, i.e. they have Tr(ρ |α, β〉 〈α ′ , β ′ |) ̸ = 0 for α ̸ = α ′ .Later in this work we will demonstrate how such modes prevent dissipative freezing.

The Monte-Carlo Unravelling of the GSKL Equation and Dissipative Freezing
In order to fully introduce dissipative freezing we must describe the Monte-Carlo unravelling of the GSKL equation in Eq. ( 1).This is an unravelling of the master equation in Eq. ( 1) into a stochastic equation of motion at the level of pure states [31][32][33].After statistical averaging of the independent solutions to this equation of motion, known as quantum trajectories, the exact dynamics of the density matrix of the system is guaranteed to be recovered.We emphasize that this unravelling of the GSKL equation is not just a mathematical formalism.If the action of the jump operator can be ascribed to a property of the environment which is measurable with a POVM (positive-operator-valued-measure), then it is possible to physically realise the stochastic equation of motion which governs the trajectories by subjecting the environment to continuous measurement1 [32,[34][35][36].As a result, the phenomenon of dissipative freezing that we will describe is relevant from both a mathematical and a physical perspective.
Given an initial pure state of the system, |ψ(0)〉, all independent solutions to the stochastic equation of motion which unravels Eq. ( 1) can be written, to first order in the propagation timestep ∆t, in the following form where i is used to index the given trajectory, N i is a normalisation constant, t = n∆t is the time with n ∈ + 0 and with n |ψ (i) ((n − 1)∆t)〉 and the operator X (i)  n which propagated it from time (n − 1)∆t to time n∆t is determined stochastically from the set S with the respective 'draw' probabilities being Figure 1: Dissipative Freezing: The Hilbert space of an open quantum system fragments into a series of subspaces in the presence of a strong symmetry operator A. For an initial state in a superposition of states in these different subspaces the long-time dynamics of a single quantum trajectory will almost always involve a 'freezing' into one of these subspaces selected at random.
The evolution of a single trajectory is thus a Markovian process: the choice of operator X (i) n is dependent only on the previous state |ψ (i) ((n − 1)∆t)〉.
The density matrix of the system at time t -i.e. the solution to Eq. ( 1) -can then be approximated as a sum over N independent trajectories: ρ(t) ≈ 1 N N i=1 |ψ (i) (t)〉 〈ψ (i) (t)| and in the limits ∆t → 0 and N → ∞ this approximation becomes exact [32].The unravelling of Eq. (1) we have just described is accurate to first order in ∆t.We opt for this first order unravelling in this paper as it is illustrative and all of our arguments are valid for arbitrarily small ∆t.Unravellings which are of a higher order in ∆t involve letting the wavefunction evolve freely under exp(−iH eff ∆t) until its norm decays below a certain threshold, at which point it is acted on with a randomly chosen jump operator.Now that we have introduced the Monte Carlo unravelling of the GSKL equation, we can explicitly describe the phenomenon of dissipative freezing: Dissipative Freezing: If an open system governed by the GSKL equation possesses a strong symmetry A then, regardless of the initial state of the system, in the limit t → ∞ a given quantum trajectory will, generically (i.e.outside of just a few known exceptional cases which will be specified in Section 2.5), only have non-zero coefficients in one symmetry subspace.
We should emphasize that the given subspace into which freezing occurs will vary with each trajectory, allowing the system to recover the full density matrix -which could have nonzero coefficients in multiple diagonal blocks -upon averaging.The probability of freezing to a given subspace will therefore be given by the initial norm of the wavefunction in that subspace.Figure 1 provides an illustration of the phenomenon of dissipative freezing.
In the case when the initial state of the system is already confined to one symmetry subspace, then the phenomenon of dissipative freezing is immediate and trivial.Each trajectory will always remain in that subspace as the operators X (i) k are block-diagonal and cannot map it out of the given subspace.This also means that once a trajectory is frozen, it will remain frozen for any future times.When the initial state of the system is in a superposition of states in different subspaces, then dissipative freezing is non-trivial.It implies a dynamical breaking of the strong symmetry at the level of trajectories -with the system randomly selecting one of the (possibly) degenerate subspaces of A in the long-time limit.As we will see, such a symmetry breaking is a direct consequence of the presence of a strong symmetry and means that the system preserves the symmetry at the ensemble level by breaking it at the level of pure states.

Recovering a Block-diagonal Steady State
The simplest case to be made for dissipative freezing can be understood from the form of the steady state.Assuming there are no inter-sector traceless non-decaying eigenmodes of the Liouvillian, then this steady state can always be written in block diagonal form just like the jump operators and Hamiltonian.We know that, if ∆t → 0 in our trajectories unravelling, upon averaging over an infinite number of quantum trajectories in the long-time limit we must exactly recover this steady state, i.e.
where the c α,β (t) are the coefficients of the ith trajectory at time t for the basis vector |α, β〉.Equations ( 8) and ( 9) therefore combine to put a constraint on the coefficients of the trajectories in the long-time limit via ensuring the steady state is block-diagonal in the long-time limit.It is clear that if dissipative freezing occurs then Eq. ( 10) is satisfied as there is only one α for which the trajectory coefficients in the long-time limit will be non-zero, implying c Without dissipative freezing, Eq. ( 10) requires the product c This is very restrictive.The symmetry subspaces are essentially independent open systems with the total Hamiltonian and jump operators being direct sums of the Hamiltonian and jump operators for each subspace.Thus, despite being independent, the subspaces constrain each other via Eq.(10).A given trajectory initialised in a superposition of subspaces must obey this constraint as t → ∞ and also ensure coefficients in the same subspace ensemble average to the appropriate steady state elements, i.e.
The more subspaces that are included in the problem (and this number is arbitrary considering they are just independent open systems), the more constraints appear on the trajectories as t → ∞.Dissipative freezing is a natural way of satisfying this constraint, independently of the number of trajectories used to approximate the steady state, the number of subspaces and the structure of the subspaces.
Moreover, dissipative freezing means we cannot use a single trajectory to simultaneously infer information (however approximate) about the steady state properties of multiple subspaces.Given their complete informational independence (if the steady state is block-diagonal it contains no information -coherences -involving both subspaces simultaneously) this makes sense.We must instead run multiple trajectories, with each providing information about a separate subspace, in order to infer such information -just as if we were to simulate these independent subspaces separately.

Long-Products of Matrices
For our second insight, and an illustration of how dissipative freezing emerges as a dynamical process, we consider the explicit time evolution of a single trajectory.To help with our analysis, we introduce the projection operator P α = d α β=1 |α, β〉 〈α, β| which projects into the subspace spanned by the eigenvectors of the strong symmetry A with eigenvalue λ α .For a given trajectory at time t we can form the projection into the α subspace as |ψ (i)  α (t)〉 = P α |ψ (i) (t)〉.Due to the block-diagonal nature of the Hamiltonian and jump operators, it follows from Eq. ( 5) that where The norm of the wavefunction in Eq. ( 11) can then be interpreted as the probability p (i) (α, t) of the trajectory i being observed in the subspace α at time t.We scale this by the normalisation constant N i (which is the same for each subspace and so irrelevant to our analysis) to define the un-normalised weight In the limit t → ∞ Eq. ( 12) is the expectation value of an infinitely long product of matrices.Long products of matrices have been studied extensively in mathematics and physics due to their relevance to dynamical systems and chaos [37][38][39].When building up repeated matrix products such as α,1 , the quantities of interest are usually the singular values σ (i)  α,q (n) := σ q (A (i) α (t)) of the product with q indexing them from largest to smallest.These singular values are intrinsically related to the powers of the matrix norm ||A (i)  α (t)|| q whose logarithm is often known as the qth Lyaponuv exponent.In general, such quantities will grow or decay exponentially with n [40,41], causing the largest singular value σ (i) α,1 (n) to become dominant when taking expectation values and allowing us to write where P(σ ) is a projector to the subspace (not to be confused with the subspaces indexed by α) spanned by the (possibly degenerate) eigenvectors of (A (i)  α (t)) † A (i) α (t) with the associated eigenvalue σ Provided the initial state in a given pair of symmetry subspaces α 1 and α 2 has a non-zero weight in their respective highest singular-value subspaces then the ratio of the probabilities p (i) (α 1 , t) and p (i) (α 2 , t) will be proportional to the ratios of the corresponding singular values, i.e.
If these singular values each vary exponentially with n, however, then the ratio of probabilities will clearly either vanish or diverge in the limit t → ∞ depending on whether the growth rate of the leading singular value is largest in α 1 or α 2 .One symmetry subspace thus becomes completely dominant and the wavefunction effectively vanishes in the other subspaces.
Freezing thus occurs between all symmetry subspaces with distinct growth rates as a natural consequence of the exponentially growing nature of A (i) α (t).We should emphasize that our statement of an exponential growth in the leading singular values σ (i)  α,q (n) is intended to be that of an average change by orders of magnitude in linear time, i.e. the quantity 2 n ln(σ (i) α,q (n)) (which can be interpreted as the average growth rate of the qth singular value) does not asymptotically tend to 0 for large n.We have been careful not to assign a specific growth rate to the singular values ∀n and instead talk about the 'average' change as we are working with matrices and only under certain conditions can one assign a single growth rate ∀n to the leading singular value of long products [37].Nonetheless, the average change of orders of magnitude of the singular values of (A (i)  α (t)) † A (i) α (t) in linear time will, via Eq.( 14), be reflected in the p (i) (α, t) and cause them to differ by increasing orders of magnitude as time increases linearly.In our numerical results we will observe this growth explicitly.

Ergodicity and Quantum Systems Subject to Repeated Measurement
For a final perspective behind dissipative freezing, we will discuss its connection between established mathematical results on the ergodicity of open systems.Specifically, work in the early 2000s successfully proved the ergodicity of an open quantum system undergoing continuous measurement in the environment when there is a single unique steady state [42,43].This ergodicity guarantees that the time-average of observables associated with a single trajectory between times 0 and τ, with τ → ∞, can be used to recover the expectation values associated with the steady state.
In this work it was also realised that, if the steady state of the system is dependent on the initial state in some manner, then this ergodic theorem cannot be proven [42].The presence of a strong symmetry in an open system guarantees the existence of multiple steady states and thus a dependence of the steady state on the initial state ρ(0) via its weights in the different subspaces.It therefore follows from Ref. [42] that an ergodic theorem cannot be established in the presence of a strong symmetry -i.e. the time-averaged dynamics of a single trajectory cannot be used to recover the ensemble average over multiple trajectories.
This lack of a standard ergodic theorem for a system with degenerate steady states implies that something else must be at play in order for the ensemble average over trajectories to recover the dynamics of the open system.This was uncovered in Ref. [44], where Kuemmerer and Maassen proved that, in this degenerate case, a single trajectory will establish ergodicity solely with respect to one of the randomly selected steady states.By observing that the existence of a strong symmetry in an open system implies multiple steady states and a series of symmetry subspaces, the work of Ref. [44] provides an ergodic understanding and proof of the phenomenon of dissipative freezing.By freezing into one of the symmetry subspaces, the trajectory of an open system can establish ergodicity and -upon time averaging -recover the steady state solely within that subspace.

The Timescale of Dissipative Freezing
With an understanding of the generality of dissipative freezing, we now tackle the question of the timescales on which dissipative freezing will occur.To aid us in our analysis, we work in the channel-state duality picture, converting operators into supervectors via and constructing the the Liouvillian superoperator which acts on these supervectors.As all of the operators involved in the construction of the Liouvillian are simultaneously block-diagonal, the Liouvillian will be too, meaning it can be expressed in the form There are thus D 2 diagonal blocks in the superoperator picture and they can be indexed by the tuples (α, α ′ ) [45].
Let us now focus on a specific block (α 1 , α 2 ) and extract the corresponding We define the following spectral gap where λ 1 , λ 2 , ..., λ d α 1 d α 2 are the eigenvalues of ( L) α 1 ,α 2 .The quantity ∆ (α 1 ,α 2 ) essentially dictates the longest timescale of the dynamics of the system in the (α 1 , α 2 ) block.If α 1 = α 2 then ∆ (α 1 ,α 2 ) is guaranteed to be 0 due to the existence of at least one steady state for that symmetry sector.If α 1 ̸ = α 2 then we refer to ∆ (α 1 ,α 2 ) as the inter-sector spectral gap and it is not necessarily zero.In fact, the assumption that the steady state of the open system is blockdiagonal is equivalent to the statement ∆ (α 1 ,α 2 ) > 0 ∀α 1 ̸ = α 2 , meaning that all inter-sector coherences will decay away with time.
Now we consider the dynamics of the density matrix of the system in this sector, which is given by P α 1 ρ(t)P α 2 .We will hereon assume t is sufficiently large such that all decay modes other than the longest-lived one have decayed away to give where λ is an arbitrary purely imaginary number, c = 〈〈ρ L ||ρ(0)〉〉 and ||ρ L 〉〉 and ||ρ R 〉〉 are the left and right eigenvectors of ( L) α 1 ,α 2 corresponding to the longest-lived mode.We can work in the trajectories picture here by substituting Eq. ( 9), giving us lim with Taking the norm of both sides we can invoke the triangle inequality to arrive at Now we define the quantity p(i) (α, t) = β |c α,β (t)| -which is larger than the standard probability p (i) (α, t) -sum over β and β ′ and use the fact β,β ′ |ρ | > 1 as the right eigenvectors are normalised to arrive at where the overbar represents statistical averaging.The product p(α 1 , t)p(α 2 , t) which emerges here is in fact a very natural quantity with which to analyse dissipative freezing as, for t → ∞, it vanishes in its presence and remains finite in its absence.Equation ( 23) dictates a bound on the rate at which it can decay and tells us that, assuming all but the longest-lived eigenmode(s) has decayed away, it cannot fall below a certain value until a certain time.For our numerics in the next section we therefore define the 'freeze-time' as the earliest time t f where p(α 1 , t)p(α 2 , t) < ε and ε is a finite, small number, ε ≪ 1.It is the time at which we consider freezing to have occurred between the two subspaces α 1 and α 2 .Following Eq. ( 23), for ε sufficiently small and N sufficiently large, this definition freeze-time is lower-bounded as providing us with an approximate timescale for observations of freezing.
Along with the bound we are able to derive on it, this definition of the freeze-time is reasonable because for a single trajectory initialised solely within the subspaces α 1 and α 2 we have that, on average, one of p (i) (α 1 , t f ) or p (i) (α 2 , t f ) will be on the order of ε 2 and thus neglibible if ε is chosen to be sufficiently small.Such a result can be seen by observing p(i) (α, t) 2 > p (i) (α, t) and p (i) (α 1 , t)+ p (i) (α 2 , t) = 1 at all times and that, at the freeze-time, We remark that it might be more desirable to achieve an inequality involving p(α 1 , t)p(α 2 , t) than p(α 1 , t)p(α 2 , t).Without making stricter assumptions about the Liouvillian, however, we are unable to achieve this and thus focus on p(α 1 , t)p(α 2 , t) when analysing the 'freezing time' in our numerics.

Exceptions to Dissipative Freezing
There are 'exceptional' cases where the conservation of 〈A〉 and the recovery of the appropriate steady state is achieved within the trajectories formalism without dissipative freezing occurring.Here we identify the two cases we are aware of and which follow naturally from the mathematical arguments of Section 2.3.

Similar Symmetry Subspaces
Consider two subspaces α 1 and α 2 which are of the same dimension and where the following is true with U being an arbitrary d α 1 × d α 2 unitary matrix and (O) α i indicating that we have extracted the d α i × d α i submatrix from O that corresponds to the α i diagonal block.Eq. ( 25) essentially states that the relevant matrices (Hamiltonian and jump operators) for the two subspaces are identical up to an arbitrary change of basis and a phase factor on the jump operators.If this is true, then freezing cannot occur between subspaces α 1 and α 2 .This is straightforwardly proven by observing that (A (i) (t) α 2 at all times and thus the weights in the two subspaces cannot diverge with respect to one another.We term the symmetry subspaces 'similar' in this case and emphasize that the steady state is still block-diagonal in this case, i.e. ∆ (α 1 ,α 2 ) ̸ = 0 ∀α 1 ̸ = α 2 and p(α 1 , t)p(α 2 , t) will decay to 0 in the long-time limit as per Eq.( 23) and Eq. ( 10) is fulfilled without dissipative freezing.
We emphasize that, although we are currently unaware of any, other forms of similarity than that expressed in Eq. ( 25) may exist between pairs of subspaces which prevents freezing.

Inter-sector Traceless Modes
The second, less trivial, case occurs when there are 'inter-sector' traceless non-decaying eigenstates of the Liouvillian, i.e.
By 'inter-sector' we mean that such states possess coherences |α 1 , β〉 〈α 2 , β ′ | between states in different symmetry subspaces and do not decay away in the long-time limit.These non-blockdiagonal modes imply, via Eq.( 23), p(α 1 , t)p(α 2 , t) cannot decay to 0 for an appropriate choice of initial state and therefore freezing cannot occur.We also note that Eq. ( 26) is equivalent to stating that Eq. ( 10) is not valid.In general, determining whether Eq. ( 26) holds requires diagonalisation of the Liouvillian and is computationally expensive -although we remark that Eq. ( 26) is guaranteed to hold if one can identify an operator, known as the strong dynamical symmetry [10,45], which satisfies [H, S] = −λS and [L j , S] = [L † j , S] = 0 where λ ∈ .As we will see in our numerical examples, however, the dynamics of a single quantum trajectory over a finite-time can provide clear evidence of Eq. (26).
Physically, the traceless modes which Eq. ( 26) implies represent information and coherences between symmetry subspaces which are protected from the system's interaction with the environment.These modes play an important role in the emergence of non-stationarity in open systems, [10], quantum synchronisation [17,18] and quantum information processing [46].
In the following section we will illustrate dissipative freezing in a number of numerical examples of open quantum systems.We will also observe the exceptions we have just described.

Numerical Details
For all of the following examples, we evolve our open system using the first-order trajectories routine described in Section 2.2.To minimise any incurred errors in the trajectory routine, we set ∆tω = 0.0025 throughout, where ω is the central energy scale in the Hamiltonian and will be specified for each example.In terms of initial states, we will opt for two different types depending on what we wish to demonstrate • I. Equal Weight: The initial state is initialised across two selected symmetry subspaces α 1 and α 2 .Its weight in each subspace is the same and, within each subspace, it is an equal superposition of all basis states Figure 2: Freezing for a 'random' Liouvillian constructed from a random blockdiagonal Hamiltonian H and a single random block-diagonal jump operator L 1 .
There are four symmetry sectors, each of dimension 4 and, in each diagonal block, Hamiltonian H consists of a random Hermitian matrix scaled by a factor of ω whilst the single jump operator consists of a random operator drawn from the complex Ginibre ensemble [47].We set γ 1 = 4ω.For plots a-b the initial state is identical and of the form of Eq. ( 28) whilst for plot c it is of the form of Eq. ( 27) with the two chosen subspaces α 1 and α 2 being arbitrary.a) Logarithm of the unnormalised weight of the wavefunction in each subspace versus time for a single trajectory.Inset) Normalised weight of the wavefunction versus time.b) Logarithm of the largest singular values of the matrix A α (t) in each subspace α for the same single trajectory.Inset) Logarithm of the ratio of the largest and second largest singular value in each subspace -values are not plotted once they can no longer be accurately determined.c) Probability density function for the individual freeze-times t f ω for 10 4 trajectories.The two different colours correspond to two independent realisations of the random Liouvillian -the initial state remains the same between them.Inset) Cumulative density function for the freeze time.We set ε = 10 −6 .
• II.Random: In each symmetry sector the initial state is set to a randomly chosen unit vector, with β α a randomly chosen integer in the range (1, d α ).
In our numerics, as per the discussion in Section 2.4, we will also be calculating the quantity p(α 1 , t)p(α 2 , t) and defining freezing as having occurred when it falls below a threshold value ε ≪ 1.We will specify the value of ε used for the given example and demonstrate the bound derived in Eq. ( 24) as well as other properties that the freeze-time obeys.

Example: Random Matrices
As our first example, we consider the Lindblad equation, as in Eq. ( 1), with Hamiltonian and Lindblad operators which are block-diagonal random matrices -i.e. the off-diagonal blocks are set to 0 and the diagonal blocks of both the Hamiltonian and Lindblad operators consist of matrices drawn from the space of random Hermitian matrices and the complex Ginibre ensemble [47] respectively.This open system is then guaranteed to possess a strong symmetry in the form of an operator which is a different scalar multiple of the identity in each block.
In Fig. 2a-b we plot the evolution of the weights (both normalised and un-normalised) of a single trajectory in each symmetry subspace, along with the corresponding leading singular value of the matrix A α (t) defined in Eq. ( 12) -dropping the i superscript here as we are focused solely on a single trajectory.We use a logarithmic y-axis in order to expose the Figure 3: Freezing for coupled spin 3/2 particles governed by the GSKL equation with the Hamiltonian defined in Eq. ( 29) and a single jump operator L = s z a .We set γ = 3ω.The index α = 1, 2, ..7 refers to subspaces of total z magnetisation −3/2, −1, ..., 1, 3/2 respectively.a) Logarithm of the un-normalised weight of a single trajectory -initialised in the form of Eq. ( 28) -versus time in each of the magnetisation subspaces.Inset) Normalised weight for tω ≤ 50.b) Probability density of the freeze-time t for 10 4 trajectories and for trajectories initialised in the form of Eq. ( 27) in the indicated subspaces with γ = 3ω.c) Average freeze-time (black dots) vs γω on a log-log scale for trajectories initialised in the form of Eq. ( 27) with α 1 = 2 and α 2 = 3 respectively.Red dotted line shows the function 40.2/∆ (α 1 ,α 2 ) whilst the black solid line shows the lower bound calculated via Eq.( 24).We set ε = 10 −6 .exponential growth of these quantities.Freezing occurs into the subspace which has the highest overall growth rate and the normalised weight in this subspace saturates to unity in the long-time limit.The exponential growth in the unnormalised weights coincides with an exponential growth of the leading singular values of the evolution matrix A α (t) in each subspace α.We also observe that the ratio of the largest and second largest singular values associated with this matrix is increasing at an exponential rate, justifying the approximation in Eq. ( 13).This growth is intrinsically tied to the initial state being propagated under an increasingly long-product of matrices.
In Fig. 2c we then plot histograms of the freeze times when measured over 10 4 trajectories for two different realisations of the random Liouvillian -where a given realisation is constructed by drawing the block-diagonal Hamiltonian and jump operator from their respective ensembles.The initial state is defined in Eq. (27).The two histograms of freeze-times display qualitatively different distributions even though the matrices are drawn from the same underlying ensemble both times.This is due to the fact that the freeze-time will be strongly controlled by the spectrum of ( L) α 1 ,α 2 , which dictates the dynamics of coherences between the α 1 and α 2 subspaces.As per Eq. ( 23), the system cannot freeze until all these off-diagonal modes have died away sufficiently.The longest-lived of these modes will thus play a crucial factor in the freezing statistics and, in this example, these modes and their eigenvalues will clearly change for each draw of the block-matrices from their ensembles.In the next example we will consider a non-random, more physical Liouvillian, allowing us to straightforwardly calculate the gap ∆ (α 1 ,α 2 ) and directly correlate it with the corresponding freeze-times and our lower bound.

Example: Coupled Spin 3/2 Particles
For this example we move onto a more physical setup consisting of two coupled spin 3/2 particles -a and b -with particle a undergoing continuous dephasing along the z axis.The local Hilbert space dimension is 4 and the Hamiltonian reads where s α a/b is the spin-3/2 operator acting in the α = x, y, z direction on either spin a or b.The dephasing is described by the single jump operator L = s z a .This system possesses a strong symmetry in the form of the total z magnetisation A = s z a +s z b and thus the Hamiltonian and jump-operator can be simultaneously block-diagonalised into seven subspaces of dimensions 1, 2, 3, 4, 3, 2, 1 and with respective total z magnetisations −3/2, −1, −1/2, 0, 1/2, 1 and 3/2.Whilst there are seven subspaces, we should not expect freezing between all of them.Specifically, the pairs of subspaces with the same absolute magnetisation are 'similar' (as per the definition in Eq. ( 25)), preventing freezing.We emphasize that are no non-decaying traceless modes of the Liouvillian here and Eq. ( 10) for these pair of subspaces even though dissipative freezing does not occur.Dynamically, this follows from the fact that, apart from the initial state, the only difference in the dynamics in the two subspaces is that when a jump occurs the coefficients of the trajectory in the negative magnetisation subspace pick up a factor of e iπ compared to the positive one.After a certain time, the total phase picked up will thus be nπ, where n is the number of jumps.The summation in the expression i c α ′′ ,β ′′ (t) can thus be split into two, one involving trajectories where n is odd and one where n is even.In the long-time limit, the probability that n is even or odd should be equal and so the two terms should be identical, aside from a scale factor of e iπ , which causes them to cancel out and the ensemble-average to equal to 0 without freezing occurring.
In Fig. 3 we plot the evolution of trajectories in this system and illustrate these results, observing freezing between subspaces of different absolute magnetisation due to the distinct exponential growth of the matrices A α (t) whilst observing an absence of freezing in subspaces with magnetisation of opposite sign due to the similarity of the subspaces.We then calculate the freezing times for trajectories initialised in a superposition of two different symmetry sectors.We observe distinctive distributions for each symmetry sector, indicating that they are specific to the form of the matrices in those sectors.We select the α 1 = 2 and α 2 = 3 sectors and analyse this further, varying γ and computing the average freezing time, the inter-sector spectral gap ∆ (α 1 ,α 2 ) associated with these sectors and the lower bound on the freezing time as per Eq.(24).
We find that our bound is obeyed for the full range of γ as expected and, whilst not being tight (which can be attributed to our application of the triangle inequality to a very large sum of complex numbers), it closely mimics the functional form of the freezing-times and demonstrates well their close relationship with the inter-sector spectral gap.For small values of γ we observe proportionality between the inverse of the inter-sector spectral gap and the freezing time.This can be understood from perturbation theory [17], with the real parts of the eigenvalues of the Liouvillian's eigenmodes all scaling proportionally (with the same proportionality constant) with γ for small γ.The inverse of these eigenvalues essentially corresponds to the timescales of the system and so the freeze-time will scale as the inverse of γ and the inverse of the spectral gap.

Example: Lossy Bosonic Chain
For our final example we consider a periodic boundary chain of L sites hosting non-interacting bosons collectively coupled to a cavity.In momentum space, the Hamiltonian reads Figure 4: Freezing dynamics for non-interacting bosons collectively coupled to a single-mode electromagnetic field.The Hamiltonian is defined in Eq. ( 30) and there is a single jump operator L = γa.We set γ = 5ω, g = 2ω, J = 2ω and the system is always initialised in an equal superposition of random states in each symmetry subspace.The index α runs over the valid momentum tuples (s 1 , s 2 , ..., s L/2 ) which specify the symmetry subspaces, the explicit relationship between α and these tuples is detailed in appendix A. We enforce a maximum photon number in our simulations of n ma x = 5.In plots a-c) we set L = 4 whilst in plots d-e) we set L = 6.a) Logarithm of the weight of a single trajectory in each of the strong symmetry subspaces verus time, initial state is of the form of Eq. ( 28).b) Logarithm of the largest singular value of A α (t) in each subspace versus time, initial state is of the form of Eq. ( 28).c) Average freeze-out time on a log-log plot for an initial state of the form of Eq. ( 27) with α 1 = 1 and α 2 = 3. Red dashed curve shows 98.3/∆ (α 1 ,α 2 ) whilst black solid line corresponds to the lower bound in Eq. ( 24).We set ε = 10 −8 .d) Logarithm of the weight of a single trajectory in each of the strong symmetry subspaces versus time, initial state is of the form of Eq. ( 28) Inset) Dynamics of the normalised weight for tω ≤ 40.e) Absolute values of the elements of the matrix p(α, t)p(α ′ , t) for tω = 250 and averaged over 10 3 trajectories, initial state is the same for all realisations and of the form of Eq. ( 28).
where a and a † , respectively annihilate and create photons in the cavity whilst b k and b † k respectively annihilate and create bosons in the chain with momenta k ∈ { 2π L , 4π L , ..., 2πL L }.We also define n k as the boson number operator for momentum k, i.e. n k = b † k b k .Photon loss is introduced into the cavity via a single jump operator L = a and we fix ourselves to halffilling, i.e. the total number of bosons is always L/2.The master equation formed from these operators was studied in Ref. [26] in the context of dissipative phase transitions of a cavityimmersed quantum gas.Importantly this setup possesses a series of L/2 commuting strong symmetry operators of the form S k = n k + n k+π with k ≤ π.
In Ref. [26] it was observed that, for the initial states considered, the individual trajectories in this system always freeze to a symmetry subspace.Here we unravel this phenomenon of dissipative freezing using the intuition we have built.Notably we also describe how, beyond a certain system size and for certain initial states, dissipative freezing is not guaranteed due to the fact that the inter-sector gap closes for certain pairs of subspaces, i.e ∆ (α 1 ,α 2 ) → 0. In Fig. 4c we again find the freezing times obey the bound we derived in Eq. ( 24).Notably as γ becomes larger this bound becomes tighter.This can be understood from the fact that for increasing γ the Hamiltonian's role in the dynamics becomes increasingly neglibible and the trajectory propagator tends towards a matrix whose elements are always positive real values.The triangle inequality used to derive our bound thus becomes tighter as γ increases and so does our lower-bound on the freeze time.It is difficult, however, to reach a numerical regime where this bound saturates as it requires using a very large value of γ (which requires a smaller ∆t to maintain accuracy), reaching exceptionally long timescales and maintaining signficant precision in the vanishingly small value of p(i) (α 1 , t)p (i) (α 2 , t).In Fig. 4c we also observe direct proportionality between the freeze-time and the inverse of the inter-sector gap for small γ -the regime in which γ acts as a rescaling of the system's timescales.
It might be natural to assume that freezing will occur in this Bosonic setup for all system sizes.In Fig. 4c-d we treat a system for L = 6 and observe that this is not the case, with freezing only occurring between certain subspaces.The origin of this lack of freezing are the traceless modes of the Liouvillian, ∆ (α 1 ,α 2 ) → 0, for L = 6 and certain pairs of subspaces.In our numerics we also observe that these modes also appear for L > 6 and independently of the maximum number of photons that we enforce.For L = 6 these modes connect the symmetry sectors where s 1 + s 2 is the same, which likely is related to the fact that the cos(k) term in the Hamiltonian only differs by a minus sign for the two separate (i.e.not simply related by a shift of +π) momentum modes k = 2π/3 and k ′ = 4π/3 which correspond to s 1 and s 2 respectively.When L ≥ 6 we can always identify such pairs of modes where cos(2πi/L) = − cos(2π j/L) and i, j ∈ 1, 2, ..., L/2.For L = 4 this is not possible.
These traceless modes are clearly manifest as non-zero off-diagonal elements in the ensemble-averaged matrix p(α 1 , t)p(α 2 , t) plotted in Fig. 4e.Following Eq. ( 23), for sufficiently long times we can be confident this matrix provides a clear picture of where the traceless non-decaying modes of the system reside (i.e.where the inter-sector spectral gap is 0).Even a single trajectory initialised across all subspaces, however, allows us to infer this information.If there are traceless modes and the freezing time diverges, then any given trajectory should have a negligible probability of freezing on any finite time-scale.
Figure 4d explicitly demonstrates this information being inferred at the single-trajectory level: all pairs of subspaces which share a traceless mode have weights of the same order of magnitude for all time, in direct agreement with the ensemble averaged picture provided by Fig. 4e.This single-trajectory picture thus immediately points to a steady state which is not block-diagonal for initial states containing coherences which overlap with these modes.The freezing time thus diverges but, as opposed to the previous setups, this happens for all finite γω and does not require it to to tend to 0 or ∞.
Importantly, our results in this paper show we do not need to wait until t → ∞ to infer this lack of freezing at the single trajectory level.Figure 4 pictures a clear manifestation of these traceless modes after a finite period of time, with the probabilities in the associated subspaces saturating to a constant, non-zero non-unity value.Given the identification of a strong symmetry and the corresponding symmetry subspaces, the ability to utilise a single trajectory to detect the existence of these traceless modes is a much more computationally efficient method than having to resort to ensemble averaging or memory-intensive calculations in the superoperator picture.
We note that these traceless modes and the resultant absence of freezing between certain subspaces was not observed in Ref. [26] as the initial states used did not span the necessary subspaces to have an overlap with these modes.

Conclusion and Outlook
In this work we have provided significant mathematical and physical insights into understanding how the recently-observed phenomenon of dissipative freezing is a general consequence of the presence of a strong symmetry in an open system.The associated symmetry-breaking in the trajectories picture is the way the system ensures symmetry is preserved at the ensemble level.We have also introduced a range of examples where dissipative freezing can be directly observed and quantified.
In these examples, we showed how the freezing statistics are dependent on the lowest lying eigenvalue in the 'off-diagonal' superoperator subspaces.When the real part of this eigenvalue vanishes, then the Liouvillian possesses traceless non-decaying modes, immediately implying an absence of dissipative freezing for any single given trajectory.We have thus identified a computationally efficient method for identifying the existence of such modes.
We envisage a number of future avenues of research stemming from this work.Firstly, identifying whether further exceptions to this phenomenon are possible or whether those identified in this work -traceless non-decaying modes and similar subspaces -are the only possible exceptions, would greatly further our understanding of the effect of symmetries on trajectories in open quantum systems.
Secondly, in this work we have adopted the approach of being as general as possible and making no assumptions about the microscopic details of the Liouvillian other than its possession of a strong symmetry.We anticipate that by making stronger assumptions about the details of the system, the bound in Eq. ( 24) can be improved and the general dependence of the the freezing properties of the system on its spectral statistics can be quantified furtherleading to novel ways to determine the spectrum of open systems from within the memory efficient trajectories formalism.
Lastly, and by no means least, an important further pursuit would be to quantify the role that symmetry-breaking perturbations play in the breakdown of dissipative freezing.Such perturbations would lead to an opening of the Liouvillian's gap in the full superoperator space and quantifying this through the freezing statistics could provide new insights into dissipative phase transitions in open systems. of Excellence ' Advanced Imaging of Matter' of the Deutsche Forschungsgemeinschaft (DFG) -EXC 2056 -project ID 390715994.We are also grateful to Berislav Buča for discussions on dissipative freezing and for suggesting the coupled model.

First, as
the strong symmetry operators all commute, there clearly exists a simultaneous block-decomposition of the Hamiltonian and jump-operator.Specifically, the corresponding subspaces or blocks can be indexed by the unique tuples (s 1 , s 2 , ..., s L/2 ) where s i corresponds to the number of bosons in momentum modes k = 2πi L and k = 2πi L + π, i.e. s i = 〈n k + n k+π 〉.If the total number of bosons is N then it follows that the number of unique tuples (or subspaces) is D = (N +L/2−1)!N !(L/2−1)![48] which grows exponentially with L for N = L/2.An initial state initialised in a superposition of states with different tuples is therefore expected to undergo the phenomenon of dissipative freezing -which is exactly what we observe in Fig. 4a-b for a small L = 4 site chain with the initial state having a non-zero weight in all subspaces.The characteristic exponential growth in the weights and singular values is observed just like in our earlier examples.