Emergence of PT-symmetry breaking in open quantum systems

The effect of PT-symmetry breaking in coupled systems with balanced gain and loss has recently attracted considerable attention and has been demonstrated in various photonic, electrical and mechanical systems in the classical regime. Here we generalize the definition of PT symmetry to finite-dimensional open quantum systems, which are described by a Markovian master equation. Specifically, we show that the invariance of this master equation under a certain symmetry transformation implies the existence of stationary states with preserved and broken parity symmetry. As the dimension of the Hilbert space grows, the transition between these two limiting phases becomes increasingly sharp and the classically expected PT-symmetry breaking transition is recovered. This quantum-to-classical correspondence allows us to establish a common theoretical framework to identify and accurately describe PT-symmetry breaking effects in a large variety of physical systems, operated both in the classical and quantum regimes.

The effect of PT -symmetry breaking in coupled systems with balanced gain and loss has recently attracted considerable attention and has been demonstrated in various photonic, electrical and mechanical systems in the classical regime. Here we generalize the definition of PT symmetry to finite-dimensional open quantum systems, which are described by a Markovian master equation. Specifically, we show that the invariance of this master equation under a certain symmetry transformation implies the existence of stationary states with preserved and broken parity symmetry. As the dimension of the Hilbert space grows, the transition between these two limiting phases becomes increasingly sharp and the classically expected PT -symmetry breaking transition is recovered. This quantum-to-classical correspondence allows us to establish a common theoretical framework to identify and accurately describe PT -symmetry breaking effects in a large variety of physical systems, operated both in the classical and quantum regimes.
The breaking of parity and time-reversal (PT ) symmetry has been widely studied in dissipative systems with an exact balance between gain and loss [1][2][3][4][5]. Owing to this symmetry, the dynamical matrix describing such systems may exhibit a purely real eigenvalue spectrum, despite a constant exchange of energy with the environment. As the dissipation rates are increased above a critical value, at least one pair of eigenvalues becomes purely imaginary and the corresponding eigenmodes no longer exhibit the symmetry of the underlying equations of motion. Over the past years, this effect has attracted considerable attention and has been demonstrated in various optical [6][7][8], electrical [9] and mechanical [10] settings. Apart from purely fundamental interest, this mechanism also has many important practical consequences, for example, for the operation of multi-mode lasers [11,12], enhanced measurements [13][14][15][16], the bandstructure of dissipative lattice systems [17] or energy transport at macroscopic [18] and microscopic [19,20] scales.
In connection with PT -symmetric systems it is common to use the terminology of non-Hermitian 'Hamilton operators'. However, the effect described above is a priori only defined for classical systems that can be modeled in terms of a complex-valued dynamical matrix [21]. Indeed, in a full master equation (ME) formulation of open quantum system [23], there is no such transition between purely real and purely imaginary eigenvalues of the corresponding Liouville operator. Also, at a microscopic level, the time-reversal equivalence between loss and gain is broken by quantum fluctuations [24][25][26][27]. Therefore, it is still an unresolved question how to formally define PTsymmetry for dissipative quantum systems [28] and if the breaking of this symmetry can exist at all at a microscopic level [27]. In several previous studies this question has been addressed by looking at coupled quantum oscillators [15,[24][25][26][27][29][30][31][32][33] or bosonic atoms [34] with gain and loss. In such systems, the symmetry-breaking effect can still be observed in the dynamics of the mean amplitudes, which simply reproduce the classical equations of motion, while quantum effects lead to increased fluctuations. However, these findings cannot be generalized to systems with a finite dimensional Hilbert space and they also provide no insigths into the stationary states of PT -symmetric quantum systems [19,26], which do not exist in purely linear models.
In this Letter we introduce a symmetry transformation for Liouville operators, which extends the conventional definition of PT symmetry to arbitrary open quantum systems. We show that under very generic conditions, the existence of this symmetry implies that the steady state of the system can be tuned between a fully symmetric and a symmetry-broken phase. While the change from one to the other limiting state is always continuous, it becomes more and more pronounced as the dimension of the Hilbert space is increased, and a sharp PT -symmetry breaking transition emerges in the semiclassical limit. This quantum-to-classical correspondence allows us to establish a unified theoretical framework for analyzing PT -symmetry breaking effects in a wide range of physical systems and to identify characteristic properties and experimentally observable features that are common to all of them.
PT -symmetric quantum systems.-We consider a generic bipartite quantum system with total Hamiltonian H. The two subsystems, A and B, have the same Hilbert space dimension, d, and are subject to dissipation described by the local jump operators c A and c B , respectively. The ME for the system density operator ρ can then be written as ( = 1) [23] where L ≡ L[H; c A , c B ] is the Liouvillian superoperator. The first term in Eq. (1) describes the evolution of a quantum state under the action of the non-Hermitian This part does not conserve the norm of the state and thus the recycling terms ∼ cρc † must be added to obtain trace-preserving dynamics. Given the decomposition of a ME in Eq. (1), it is tempting to define PT -symmetric quantum systems in analogy to the classical case [4,5], namely as open quantum systems where (PT )H eff (PT ) −1 = H eff . Here P is the parity operator with P(A ⊗ B)P −1 = B ⊗ A and T iT −1 = −i. However, H eff has only negative imaginary parts because the norm of a state evolving under H eff always decreases. Thus, this symmetry relation can only be satisfied in closed systems. The same is then also true for the eigenvalues of the full Liouville operator L whose real part must always be negative or zero. In conclusion, while there is a natural way to introduce non-Hermitian Hamiltonians in open quantum systems and even probe them via conditional measurements [35][36][37][38][39], there are no PT -symmetric (super-)operators in the conventional sense.
To extend the concept of PT symmetry into the quantum regime, it is important to keep in mind that the relevant physical effect of the T -operator is to exchange loss and gain and not necessarily to implement a full time-reversal transformation. In the simplest example of a quantum harmonic oscillator the effect of loss with rate Γ is modeled by a jump operator c = √ Γa, where a is the annihilation operator. In turn, the effect of gain with the same rate can be described by modifying the jump operator to be c = √ Γa † . Therefore, in this case we find that the transformation between loss and gain is implemented in the ME formalism by replacing the jump operator by its adjoint operator, c → c † . Guided by this explicit example, we introduce the following anti-unitary transformation for operators O, and define an open quantum system to be PT -symmetric, if the corresponding Liouvillian satisfies This condition implies that the Hamiltonian H is paritysymmetric and that the local jump operators are of the form where O can be an arbitrary dimensionless operator. We remark that this definition differs from the PTsymmetric Liouville operators introduced in Ref. [28], where, to our knowledge, the considered transformations have no immediate physical interpretation or classical correspondence. While the systems studied in In (c) and (d) we plot the corresponding dependence of the symmetry parameter ∆ defined in Eq. (7) on the ratio Γ/g. In (c) the line for S = ∞ is obtained from a Holstein-Primakoff approximation [44]. See text for more details.
Refs. [28,40,41] satisfy Eq. (4) with a redefinition of P, none of the examples discussed below exhibits the symmetry considered in these references when d > 2.
Phenomenology.-Before we return to a more general discussion of Eq. (4), let us illustrate its physical implications in terms of two simple examples: (i) Two coupled spin S = (d − 1)/2 systems with O = S + , where S + = S x + iS y is the spin raising operator, and (ii) two coupled harmonic oscillators with O = a † . In the second example we introduce a finite cutoff occupation number, i.e., a † |n = d − 1 = 0. This cutoff mimics the effect of saturation in realistic systems [19] and allows us to vary the Hilbert space dimension. In both examples we consider a Hamiltonian of the form where This Hamiltonian describes the coherent exchange of energy between the two subsystems with a strength g. The resulting Liouvil- , then satisfies Eq. (4). We calculate the steady state, ρ 0 , satisfying Lρ 0 = 0, for different ratios Γ/g and show in Fig. 1 the symmetry parameter [26] This is an experimentally observable quantity, only requiring measurements of local operators, which provides a measure for the symmetry of the system, i.e., ∆ = 0 for a parity-symmetric density operator, PρP −1 = ρ. For the current examples, ∆ represents the normalized population imbalance between the two subsystems. For small dimensions d, this parameter changes gradually from 0 to 1 with increasing Γ. This smooth variation is expected since observables of finite dimensional quantum systems cannot exhibit any non-analytic behavior. However, as the system size increases, ∆ vanishes for Γ/g < 1 in the limit d → ∞, while it retains a finite value for Γ/g > 1.
In both examples, the critical ratio is Γ/g = 1, which corresponds to the dynamical PT -symmetry breaking point of an equivalent linear oscillator system with gain and loss [4,5]. We thus conclude that PT -symmetry breaking exists even for non-harmonic and finite dimensional quantum systems, but only as an emergent phenomenon in the semiclassical limit.
To obtain better insights into the nature of the two phases, we plot in Fig. 2(a) the purity, P = Tr{ρ 2 0 }, for the steady state of the spin system. This quantity again exhibits a sharp transition around Γ = g and shows that the symmetric and symmetry-broken phases are characterized by a highly mixed and an almost pure steady state, respectively. More precisely, the scaling P (Γ → 0) d −2 implies that in the symmetric phase the steady state is close to the maximally mixed state, ρ 0 (Γ g) 1/d 2 . This indicates that for Γ < g the gain and loss processes cancel out on average while quantum fluctuations still occur with rate Γ and completely randomize the system's long-time dynamics [19,26]. In contrast, for Γ > g, the incoherent processes dominate and pump the spins into the polarised pure state, Closer to the transition point, the coherent coupling creates excitations ∼ O † A O B |ψ 0 on top of this state, which are strongly correlated. As shown in Fig. 2(b), this results in a characteristic peak in the entanglement between the two subsystems around the transition point, which vanishes again in the symmetric phase due to fluctuations. Consistent with similar features observed in saturable oscillator systems [19], this peak in the entanglement shows that even for d 1 the PT -symmetry breaking transition retains genuine quantum mechanical properties.
Existence of a fully symmetric steady state.-We will now show that the properties discussed above for specific examples are indeed a general consequence of the symmetry relation in Eq. (4). Firstly, we demonstrate that for any Liouvillian that satisfies this condition and where the spectrum of H is non-degenerate the fully mixed state, is a stationary state of L in the limit Γ → 0. To do so is the coherent evolution, and write the density operator as ρ = n,m ρ n,m |E n E m |, where H|E n = E n |E n . For Γ = 0 any diagonal state with ρ n,m = 0 for n = m is a stationary solution of the ME, but the populations p n = ρ n,n are not uniquely determined. To show that for small but finite Γ only the fully mixed state is dynamically stable, we write p n = 1/d 2 + δp n . Up to first order in Γ we then obtain [44].
We now make use of the relation Pc B P −1 = c † A , which follows from Eq. (5), and the fact that the eigenstates of H are also eigenstates of the parity operator, i.e., P|E n = ±|E n . This allows us to A ]|E n and δṗ n = 0. This result shows that for PT -symmetric quantum systems the fully mixed state is stationary in the presence of a small amount of dissipation, even when each individual jump operator c A,B would drive the system into a polarized state. Note that the analysis presented here assumed a non-degenerate spectrum, i.e., the absence of any additional symmetries S other than parity. In the general case the same arguments still hold as long as all eigenstates with different parity are separated by a finite energy gap, or, more formally, as long as [S, P] = 0 [44].
Symmetry-breaking transition.-While the existence of a fully symmetric steady state follows directly from Eq. (4), there are many trivial cases where this is also the only stationary state, for example, when O is Hermitian. Therefore, we are interested in systems where there is a competing asymmetric phase in the limit Γ → ∞. To ensure that such a phase exists we now restrict ourselves to a Hamiltonian as given in Eq. (6) and a non-Hermitian jump operator of rank d−1 with Tr{O} = 0. This implies that there are dark states |D and |D * , which satisfy O|D = 0 and O † |D * = 0. Under these assumptions we obtain the symmetry-broken phase which is fully asymmetric, ∆ = 1, and has maximal purity, P = 1. Note, however, that for observing symmetrybreaking effects it is not essential that ρ 0 (Γ → ∞) is a pure state and in [44] we discuss an examples, where the symmetry-broken state is mixed. Given the two distinct limiting phases, the remaining question is, if there is a sharp phase transition between them at a critical intermediate value Γ c . For the spin system discussed above this question can be rigorously answered in the limit S 1 by examining the stability of linear fluctuations on top of the fully polarized state. This can be done using a Holstein-Primakoff approximation [1], where the spin operators are replaced by a pair of bosonic operators, From the analytic solution [44] of this linearized model we find that the fluctuations a † a and b † b diverge at the point Γ c = g. This analysis also predicts very accurately the observed dependence of the purity and the maximum of the entanglement [2], N (Γ = Γ c ) = 1/2, as indicated in Fig. 2(a)-(b).
In general, such an analytic treatment is not possible and, in many situations, PT -symmetry breaking can occur as a smooth crossover, rather than a sharp phase transition. Nevertheless, it turns out that the appearance of a sharp transition in the limit of large d does not require any specific fine tuning of the dissipation mechanism. This point is illustrated in Fig. 2(c)-(d), where we consider a set of PT -symmetric quantum systems with randomly generated jump operators O [44]. For each individual instance, we observe the characteristic transition between the fully mixed and pure states and the asymmetric entanglement peak. These features sharpen as the Hilbert space dimension is increased. Therefore, this study demonstrates that sharp PT -symmetry breaking transitions are not restricted to simple systems with a direct classical counterpart and are expected to occur in a large range of systems that obey Eq. (4).
Discussion and conclusions.-In summary, we have introduced the symmetry relation for Liouville operators, which extends the notion of PT symmetry to bipartite open quantum systems. This definition is consistent with previous examples of linear PT -symmetric quantum systems for which the conventional definition of PT symmetry is recovered in the limit of large oscillation amplitudes. At the same time the map PT in Eq. (3) is completely general and can also be used to study PT symmetry in highly nonlinear systems or for dissipation processes that have no direct classical counterpart. The symmetry can be immediately generalized to systems with multiple jump operators, or to systems where . While this model contains only loss processes and the occupation numbers S + A S − A = S + B S − B remain symmetric for all ratios of Γ/g, the general arguments presented in this paper still hold and indeed this model exhibits a sharp transition from a fully mixed to a highly pure state [44]. The symmetry relation in Eq. (4) is thus a powerful tool to identify PT -symmetry breaking effects, even in systems where our naive intuition fails.
In this paper we have mainly focused on the steady state ρ 0 , which is determined for all parameters by the zero eigenvector of L. Traditionally, PT -symmetry breaking is discussed as a qualitative change in the system dynamics, which is associated with the appearance of exceptional points in the eigenspectrum of the dynamical matrix. This has motivated similar studies of the spectra of Liouville operators, where the appearance of exceptional points [48,49] or additional symmetries in the complex eigenvalue structure [28,40,41] have been discussed. Given the generalized definition of PT symmetry introduced here, we consider the example of two spin S = 4 systems from above, for which we plot the full Liouville spectrum below and above the transition point in Fig. 3(a)-(b). While for the two cases we don't observe any significant differences in the overall eigenvalue structure, the evolution of the observables S z A,B shown in Fig. 3(c)-(d) still undergoes the expected transition from an oscillatory to an overdamped behavior. This final example confirms our previous conclusion, namely that PT -symmetry breaking is an emergent phenomenon in the dynamics and stationary expectation values of macroscopic observables, which, in general, depend little on individual eigenvalues of H eff or L. Based on the symmetry in Eq. (4), this effect can now be studied more systematically in a large variety of physical systems to make accurate predictions for upcoming experiments in this field. Our discussion also shows that there are still many interesting conceptual questions to address. This concerns, in particular, the existence and the nature of the PT -symmetry breaking transition in higher dimensional lattices and interacting many-body system, for which no exact numerical solutions exist.

EXISTENCE OF A FULLY MIXED SYMMETRIC PHASE
In this section we detail and extend the proof for the linear stability of the fully mixed symmetric phase in the limit Γ → 0 discussed in the main text. As a starting point we write the density operator as where |E n are the energy eigenstates of H, i.e. H|E n = E n |E n . Since, according to the PT -criterion in Eq. (4) of the main text, [H, P] = 0, we may simultaneously diagonalise the parity operator and hence P|E n = ζ n |E n , where |ζ n | 2 = 1 without loss of generality.
For Γ = 0 the fully mixed state, ρ = 1/d 2 , is a stationary solution of the MEρ = L g ρ = −i[H, ρ], but this is also true for any other diagonal state. Therefore, we make the ansatz ρ n,m = δ n,m /d 2 + δρ n,m and evaluate the evolution of δρ n,m up to first order in Γ [note that We first assume that E n = E m . In this case the elements ρ n,m represent coherences between non-degenerate eigenstates and we obtain Therefore, to lowest order in Γ all these off-diagonal elements of the density matrix remain bounded and |δρ n,m | → 0 for Γ → 0. For all other matrix elements with E n = E m the coherent evolution vanishes and δρ n,m = 1 This results in a linear growth in time, unless the matrix element on the right-hand side is zero. We now make use of the relation which follows from the PT -symmetry relation for the Liouville operator. Based on this transformation we obtain and the evolution equation from above can be written as In the case of a Hamiltonian H with a non-degenerate spectrum, Eq. (S7) only applies to the populations p n = ρ n,n , in which case |ζ n | 2 = 1 and the right hand side vanishes. This is the result given in the main text.
A bit more care must be taken for Hamiltonians with degeneracies imposed by extra symmetries beyond that generated by P. Even though the populations in a given basis still remain fixed, the build-up of coherences between degenerate levels leads to a deviation from the fully mixed state. If the Hamiltonian has a symmetry, S, such that [H, S] = 0, then the states generated by applying S to |E n are degenerate. From Eq. (S7) we see that this leads to a non-identity steady state when two states |E n and |E m with the same energy have a different parity, ζ n = ζ m . However, if [P, S] = 0 then it is straightforward to see that ζ n = ζ m . Therefore, for the existence of a fully mixed symmetric phase it is in general not enough that [H, P] = 0. In addition, we require that all other non-trivial symmetries of the Hamiltonian also commute with the parity operator, at least within each degenerate subspace.
A simple example where such non-trivial symmetries play a role is the spin model described by the Hamiltonian and the PT -symmetric MEρ This model has a symmetry generated by S = S z A − S z B which does not commute with P and indeed one can show that the steady state for this model has spin-A pointing down and spin-B pointing up independent of the value of Γ/g.

LARGE Γ/g LIMIT
The Holstein-Primakoff transformation [S1] provides an exact mapping from the collective spin operators S ± , S z onto a bosonic mode with annihilation operator c. This approach allows us to find a description of the spin system discussed in the main text for the PT -broken phase valid in the large S limit. For example the spin on the B site, which undergoes decay, is close to the spin down state with S z B = −S and hence we use this as the ground state of the oscillator and obtain (S10) In the large spin limit, S → ∞, this transformation can be used to find a linear representation of the spin operators.
In the above case this is accurate in the limit c † B c B /(2S) 1 leading to the linear approximation Equivalently, for the spin in site A, which is pumped close to the spin up state, where S z A ≈ S, we find After making this transformation the PT -symmetric ME readṡ with the Hamiltonian This ME is quadratic in the bosonic operators and the resulting set of equations of motion for the variances c † i c i , c 2 i , c † i c j , etc. can be solved analytically. From the solutions for c † i c i we obtain the following magnetization of the spins in steady state .
(S16) (S17) We can clearly see that this solution breaks down for Γ ≤ g, giving the location of the PT -symmetry breaking phase transition in this model as Γ c = g.

Purity and Entanglement
For Gaussian states we can calculate the purity and entanglement negativity from the covariance matrix [S2]. Since within the Holstein-Primakoff approximation the steady-state is Gaussian we may examine these quantities to understand more about the nature of the phase. The covariance matrix for the dimer is defined as where the X i are the quadrature operators defined as . The covariance matrix has the following structure where V A contains correlations within the first site, V B those in the second site and V C the cross-correlations.
In the symmetry-broken phase described by the Holstein-Primakoff approximation the two mode squeezing in the Hamiltonian, Eq. (S14), generates entanglement between the two spins. This can be quantified by the entanglement negativity N , which reads The purity can be obtained as These are the results which are shown as the S → ∞ curves in Fig. 2 of the main text.

RANDOM JUMP OPERATORS
In Fig. 2 (c) and (d) of the main text we calculate the steady state of random PT -symmetric finite dimensional quantum systems. Here we describe how these random models are constructed. For simplicity we keep the relationship between jump operators and the Hamiltonian as described in the main text. We also wish to ensure that the jump operators have a single dark state, such that in the limit Γ → ∞ the purity P → 1.
The procedure we use is then as follows: We first create a random matrix R from the Gaussian orthogonal ensemble (GOE), i.e., a symmetric matrix with real entries which follow a Gaussian distribution [S3]. This matrix is then shifted by its lowest eigenvalue such that R = R − λ 0 I is positive semidefinite with a guaranteed zero eigenvalue. To obtain the jump operator O we then perform a Cholesky decomposition on the resulting matrix, such that O is a lower triangular matrix. Since the Cholesky decomposition for positive semi-definite matrices is not unique, we implement this step by first diagonalizing the random matrix R , with U a unitary matrix and D = diag(0, λ 1 , . . . , λ d−1 ), a diagonal matrix where λ n are non-zero eigenvalues. The diagonal matrix D can be decomposed as D = LL † , where only the first super diagonal of L † is non-zero with ( √ λ 1 , √ λ 2 , . . . , λ d−1 ). As a result the jump operator is This procedure of constructing a random jump operator ensures that most of the resulting decay rates are O(1), due to the fact that the spacing between the eigenvalues of R will follow a Wigner surmise distribution P (∆E) ∼ ∆E exp(−A∆E 2 ) [S3], meaning that there are very few almost degenerate states. By enforcing L † to only have non-vanishing elements in the first upper diagonal ensures that it is possible to observe the PT -symmetry breaking transition. This is not guaranteed in general. For example, by decomposing the diagonal matrix D in Eq. (S23) in terms of two diagonal matrices D = √ D √ D, the resulting jump operator would be Hermitian and there would be no phase transition since the trivial identity state is always a steady state of such a model.

PT -SYMMETRY BREAKING WITH MULTIPLE DISSIPATION CHANNELS
In the main text we focused on cases where each subsystem is only affected by a single dissipation channel with a unique dark state. This implies that the steady state for Γ → 0 is pure. However, this is not a crucial restriction and the definition of PT -symmetric Liouville operators can be straightforwardly generalized to multiple dissipation channels. As an illustrative example, we consider a spin model very closely related to that presented in the main text. The Hamiltonian is again but in the ME we now include two dissipation terms for each spin, Here the different jump operators are This model represents two coupled spins, where one is coupled to a positive temperature reservoir while the other is coupled to a negative temperature reservoir. For p = 1, we recover the model in the main text. In Fig. S1 we plot the symmetry parameter, the purity and the entanglement negativity for p = 0.8. In this case the symmetry-broken phase in the limit Γ → ∞ is also a mixed state. However, all the signatures of the PT -symmetry breaking transition are still clearly visible.

EXTRA UNITARY SYMMETRIES
In the main text we mention the possibility to complement the parity transformation by an addition unitary transformations, P → PU . This means the PT map defined in Eq.
satisfies the symmetry relation in Eq. (4) of the main text. As shown in Fig. 2(a), in this case the steady state of the system is always symmetric and the jump operators for the two subsystems both induce loss. Nevertheless, since [H, PU ] = 0 all the arguments presented in Sec. I above still hold and therefore we expect also in this system a transition between a fully mixed state and a polarized phase. This is confirmed by the plot of the purity in Fig. 2(b).
Indeed, for the current example this behavior can be understood from the fact that the present system differs from the one discussed in the main text only by a flip of the first spin, i.e., S ± A → S ∓ A and S z A → −S z A . Therefore the two systems are unitarily equivalent.