PT-symmetry breaking in the steady state of microscopic gain-loss systems

The phenomenon of PT (parity- and time-reversal) symmetry breaking is conventionally associated with a change in the complex mode spectrum of a non-Hermitian system that marks a transition from a purely oscillatory to an exponentially amplified dynamical regime. In this work we describe a new type of PT-symmetry breaking, which occurs in the steady-state energy distribution of open systems with balanced gain and loss. In particular, we show that the combination of nonlinear saturation effects and the presence of thermal or quantum noise in actual experiments results in unexpected behavior that differs significantly from the usual dynamical picture. We observe additional phases with preserved or `weakly' broken PT symmetry, and an unconventional transition from a high-noise thermal state to a low-amplitude lasing state with broken symmetry and strongly reduced fluctuations. We illustrate these effects here for the specific example of coupled mechanical resonators with optically-induced loss and gain, but the described mechanisms will be essential for a general understanding of the steady-state properties of actual PT-symmetric systems operated at low amplitudes or close to the quantum regime.

state energy distribution of open systems with balanced gain and loss. In particular, we show that the combination of nonlinear saturation effects and the presence of thermal or quantum noise in actual experiments results in unexpected behavior that differs significantly from the usual dynamical picture. We observe additional phases with preserved or 'weakly' broken PT symmetry, and an unconventional transition from a high-noise thermal state to a low-amplitude lasing state with broken symmetry and strongly reduced fluctuations. We illustrate these effects here for the specific example of coupled mechanical resonators with optically-induced loss and gain, but the described mechanisms will be essential for a general understanding of the steady-state properties of actual PT -symmetric systems operated at low amplitudes or close to the quantum regime. arXiv:1508.00594v3 [quant-ph] 2 Aug 2016

I. INTRODUCTION
In 1998 Bender and Boettcher [1] described a class of non-Hermitian 'Hamiltonians' that exhibit a purely real energy spectrum, a surprising fact which they attributed to the underlying PT (parity-and time-reversal) symmetry. Their observation triggered considerable interest in discrete and continuous systems with PT symmetry along with alternative non-Hermitian formulations of quantum theory [2]. Such fundamental considerations remain speculative, but there exist many classical systems in which PT -symmetric dynamics can be obtained with appropriately engineered loss and gain. This was first pointed out in the context of photonic waveguides [3][4][5][6], lattices, and resonators [7][8][9][10]. Other examples include cold atoms [11][12][13] and optomechanical devices [14][15][16].
Of particular interest in such systems is the breaking of PT symmetry, i.e., when by tuning a parameter the energy spectrum becomes complex and the eigenvectors no longer exhibit the underlying PT symmetry. This phenomenon was first experimentally observed in optical waveguides [5,6], and is currently the subject of intense experimental and theoretical research. importantly, we identify an unconventional transition from a high-noise balanced energy distribution to a parity-broken lasing state with strongly reduced fluctuations. This transition generalizes the phenomenon of PT -symmetry breaking-hitherto defined only for eigenstates-to steady-state distributions of noisy systems. The mechanisms described here occur in systems of two coupled modes as well as in multi-resonator arrays, and will thus be of relevance for a large range of PT symmetric systems operated at low amplitudes and close to the quantum regime. and |y lead to a net absorption or emission of phonons of frequency ω m (see reference [21] for more details).

II. MODEL
To motivate the following analysis by a concrete physical system, we consider a setup of two micromechanical resonators as shown in figure 1 a). The main findings discussed below, however, are more general and can be studied in other equivalent realizations, e.g., with coupled optical or microwave resonators. The mechanical resonators have a bare vibrational frequency ω m and they are mutually coupled, e.g., mechanically via the support, with strength g. In addition, optical or electrical cooling [17][18][19][20][21] and pumping [21][22][23][24][25][26][27][28] schemes are used to induce mechanical loss for one resonator and an equal amount of mechanical gain for the other. In a frame rotating with ω m , the semiclassical dynamics of the system is then described by the Itô stochastic differential equation [29,30]  α β where α and β are the dimensionless amplitudes of the pumped and cooled mode respectively (more details on the derivation of equation (1) are given in A).
The optically-induced gain and loss rates considered here are of the form where Γ is the maximal rate and √ n 0 is the saturation amplitude. The value of ν characterizes the underlying heating or cooling mechanism and will be treated here as an adjustable parameter. For the three-level scheme depicted in figure 1 b) this parameter takes the value ν = 2 [21]. Instead, for conventional laser amplification with inverted two-level systems we would obtain ν = 1 [29].
Finally, γ is the bare mechanical damping rate. Since we are interested in the PT -symmetric limit (defined below), we assume γ/Γ → 0. However, in all our calculations we retain a finite γ > 0, which describes the actual physical situation and results in a well-defined steady state for all parameter regimes.
In equation (1) the (complex) stochastic forces F ± (t) represent two independent white-noise For resonators coupled to a reservoir of temperature T the diffusion rates are D + (α) = D q (α) + 2γN th and D − = 2γN th , where N th = (e ωm/k B T − 1) −1 .
The contribution D q (α → 0) = 2Γ for the gain mode represents the intrinsic quantum noise associated with any amplification process and suggests that noise is a fundamental property of PT -symmetric systems [31][32][33].

III. PT -SYMMETRY BREAKING
By ignoring the effect of noise and assuming |α| 2 , |β| 2 n 0 , equation (1) can be written as ∂ t ψ = −iHψ, with a state vector ψ = (α, β) T and a non-Hermitian 'Hamiltonian' where sin(θ) = Γ/g. We see that despite being non-Hermitian, for Γ ≤ g both eigenvalues λ ± are real, indicating purely coherent oscillations between the two modes. In this regime of unbroken PT symmetry the eigenvectors are eigenstates of the symmetry operator, i.e., PT ψ ± = ψ ± . For Γ > g, both eigenvalues become imaginary and correspond to a gain and a loss eigenmode. In this regime θ is complex and ψ ± no longer possess the same symmetry as H. One thus speaks of a PT -symmetry-breaking transition occurring at the (exceptional) point Γ = g [2].
This conventionally studied PT -symmetry-breaking effect captures the change in the system's transient dynamics, as it is observed, e.g., in the propagation of coupled optical modes through an active medium with loss and gain [5,6]. However, a simple eigenvalue analysis does not afford any conclusions concerning its steady state. Firstly, due to an exponentially amplified mode in the symmetry-broken phase, nonlinear effects must be taken into account in order to determine the system's long-time behavior [4,[35][36][37][38]. Secondly, due to non-decaying oscillations in the PTsymmetric phase, any source of noise would heat up the system indefinitely, and within a linearized and |β ss | 2 as a function of Γ/g, and for the two relevant cases ν = 2 and ν = 1. Firstly, we observe in both plots the expected transition at Γ/g| I→II = 1, below which (phase I) the system dynamics is oscillatory with a small overall damping to |α ss | 2 = |β ss | 2 = 0 due to a finite γ > 0.
Above this transition point (phase II) the linearized system dynamics becomes unstable and both resonators reach a finite steady-state occupation number |α ss | 2 /n 0 = |β ss | 2 /n 0 = (Γ/g) 1/ν − 1 + O(γ), determined by the saturation of Γ ± (α). However, this steady state is still an eigenstate of the symmetry operator, PT ψ ss ∝ ψ ss , and contrary to our naïve expectation the system remains PT -symmetric beyond the conventional transition point.
The existence of a PT -symmetric steady-state with non-vanishing amplitude is permitted in models like the present one, where PT -symmetry is retained even in the nonlinear regime. This means that a steady state ψ ss with |α ss | = |β ss | = 0 results in an equal suppression of both the gain and the loss rate, Γ + (α ss ) = −Γ − (β ss ) [see equation (2)]. This implies that there exists a symmetric state which satisfiesψ ss = 0 for all values of Γ/g. However, as one can see in figure 2, at larger values of Γ the system eventually switches to a different, symmetry-broken state. As we discuss now, the details of this PT -symmetry-breaking mechanism and the resulting state depend on the actual form of the saturable gain, which is here determined by the parameter ν.
Considering first the case ν = 2, we find a second transition at Γ/g = 4, beyond which (phase and the steady state breaks parity, i.e., |α ss | = |β ss |. The PT -symmetry-breaking point can be obtained from a linear stability analysis for phase II and for the present model we find This result would imply the absence of symmetry breaking for ν = 1. Instead, for this case we observe a Hopf bifurcation around Γ/g ≈ 5.2, beyond which the system approaches a limit cycle and undergoes periodic oscillations with the characteristic frequency ω osc ≈ 2 g 3 (Γ − g)/Γ. The limit cycle is formed symmetrically around a central point ψ c with |α c | = |β c |, meaning that PTsymmetry is preserved on average, but not at each point in time. To distinguish this behavior from the previous case we say that in this phase (IIIw) the PT symmetry is only 'weakly' broken. In contrast to full symmetry breaking, this transition is related to the finite asymmetry induced by γ > 0, and would be hidden in the idealized models where γ = 0. The value for the transition point is independent of the precise value of γ Γ. The general expressions in equations (5) and (6) show that for all intermediate cases, 1 < ν < 2, we obtain 1 < (Γ/g) II→IIIw < (Γ/g) II→III < ∞ and PT -symmetry breaking occurs in three steps I → II → IIIw → III. More details on the derivation of the above results and the stationary phases for general ν are given in B.
While similar nonlinear phenomena are in general expected for gain-loss systems [16,37,39,40], our specific interest here is to understand the role of dynamical instabilities in the breaking of a steady-state symmetry. In particular, the above analysis shows that PT -symmetry breaking occurs even in systems where the symmetry is fully retained in the nonlinear regime and a symmetric steady-state would be permitted in principle for all parameters. Plot of the PT -symmetry parameter ∆ defined in equation (8). c) Steady-state distribution of α (red dots) and β (blue dots) in the thermal (Γ/g = 2, left plot) and in the symmetry-broken (Γ/g = 10, right plot) phases.

V. PT -SYMMETRY BREAKING IN NOISY SYSTEMS
We now show how the above picture changes in the presence of noise. For clarity we first restrict ourselves to a system which is dominated by thermal diffusion, i.e. D ± = D = 2γN th . For Γ = 0 thermal noise induces additional amplitude fluctuations of about (∆α) 2 ≈ D/(2γ) = N th , and we expect that as long as N th < n 0 the characteristic features shown in figure 2, which scale with the saturation number n 0 , will only be smeared out, but not change significantly. This is confirmed numerically (not shown) and means that figure 2 is a good representation of the steady states of the system in the weakly nonlinear or low-noise regime. Therefore, we will now address the opposite regime N th n 0 .  (1), from which we obtain the steady-state distribution P ss (α, β) for ν = 2 and N th /n 0 = 10. In the following we write α = re iθα and β = ze iθ β and make use of the fact that the system dynamics is invariant under a global phase rotation. The exact marginal distributions P ss (α) and P ss (β) are then fitted by approximate distributions of the form which allows us to extract a radial shift r 0 and z 0 and the range of fluctuations (∆α) 2 and (∆β) 2 for both modes (see C). From these values plotted in figure 3 a) we see that the thermal noise now completely washes out the features associated with the PT -symmetric phases I and II, and for a large range of Γ the system reaches a steady state (phase T), which is to a good approximation thermal, i.e., r 0 = z 0 = 0 and (∆α) 2 (∆β) 2 N th . Only after a critical value of Γ/g| T→III ≈ 7.5 are the fluctuations suddenly strongly suppressed. In this regime the system relaxes into an asymmetric coherent state with r 0 > z 0 approximately given by the amplitudes |α ss | and |β ss | given in equation (4) and (∆α) 2 , (∆β) 2 ∼ γN th Γ/g 2 , γN th /Γ 1.
Before we proceed, let us connect this transition to the phenomenon of PT -symmetry breakinghitherto defined only for individual states. To do so we introduce the PT -symmetry parameter which vanishes for a random set of states for the amplitude of the loss mode, z = |β| [41], where η z (t)η z (t ) = δ(t − t ) and (for γ → 0) The function U (z) is an effective potential for z, which is sketched in Fig. 4 a). This potential has a local minimum at z 0 = |β ss | (corresponding to the steady state given in equation (4) for N th → 0), which is separated by a finite barrier ∆U from the unstable region z > z max . In the presence of noise, a system initially located at z ≈ z 0 can escape over this barrier via thermally activated processes with a characteristic rate This rate increases as Γ is reduced and once R esc exceeds the bare damping γ, any configuration with fixed α and β is rapidly destabilized and the transition to a quasi-thermal state with strongly fluctuating amplitudes occurs. In figure 4 b) we compare the transition point Γ/g| T→III obtained from the condition γ = R esc with the numerically evaluated values for various N th /n 0 1. The plot shows that PT -symmetry breaking in the large-noise regime is very well described by this thermal activation model.

VI. QUANTUM NOISE LIMIT
As it has been pointed out above, the implementation of any gain mechanism is necessarily accompanied by a minimum amount of quantum noise, D q , which ensures, e.g., the preservation of Heisenberg's uncertainty relations in a quantum mechanical amplification process. For PT - and D q (α) as defined in equation (11). The other parameters for this plot are n 0 = 10 and γ/g = 10 −3 .
symmetric systems this implies that noise is not only an intrinsic property, but also that the symmetry in the gain and loss modes is broken on a fundamental level.
To illustrate the effect of pure quantum noise on the system's steady state we set N th = 0 and consider the quantum diffusion term only. The assumed form of D q (α) applies for the three-level gain mechanism shown in figure 1, i.e., for the case ν = 2 (see A for additional details). Similar to the thermal noise limit discussed above we represent the resulting steady-state distributions for the loss and gain modes by the distribution maxima r 0 , z 0 and the range of fluctuations ∆α and ∆β respectively. The results are shown in figure 5 for a saturation parameter n 0 = 10. We see that while the steady-state distributions for low Γ do no longer represent thermal states, there still exists a sharp transition between a high-noise PT -symmetric phase and a PT -symmetry-broken phase with strongly reduced fluctuations. Thus, we conclude that the PT -symmetry-breaking mechanism identified above is not very sensitive to the details of the noise process and will exist in both thermal and quantum noise limited systems.
Note that the semi-classical descriptions used in the current analysis is only valid for n 0 1 while for n 0 ∼ 1 where quantum effects are most important a full quantum mechanical treatment must be employed. Such an analysis is, however, beyond the scope of the current work and will be presented elsewhere.

VII. ARRAYS
Finally, to show that PT -symmetry breaking in the steady state exists also for extended systems, we generalize the analysis above to coupled resonator arrays with alternating gain and loss [42,43], as depicted in figure 6 a). In this case the amplitudes α n and β n for each unit cell obey where g is the coupling between the unit cells and F n,± (t) are independent thermal noise processes. Figure 6 summarizes the numerical results for the steady state of an array of N = 12 resonators with periodic boundary conditions and ν = 2. The observed features can be understood from the plane wave ansatz [42,44] α n = A k e ikn , β n = B k e ikn , where k = 4πj/N . This ansatz maps equation (12) onto a two-mode problem for A k and B k , which is equivalent to equation (1), but with the replacement g → g k = |g + g e ik |. This implies that the largest ratio Γ/|g k | is always achieved for the k = π mode [43], which therefore determines the symmetry-breaking properties of the array.
For the special case g = g , i.e., g k = 0, the gain and the loss modes completely decouple, and PT -symmetry breaking already occurs for Γ → 0 [9]. For all intermediate parameters the phase boundaries in figure 6 are obtained from the analytic results for the two-mode problem, but with g replaced by g − g . We see that the single-mode ansatz captures well the relevant physics both in the low-and high-noise regime. Note that, however, in the 'thermal' phase the behavior of the array can actually be much more complicated, since the system may undergo noise-induced jumps between multiple metastable configurations.

VIII. CONCLUSIONS
We In this section we outline the main steps in the derivation of the stochastic differential equation (1). To do so we first consider an interaction between a single mechanical mode and a second auxiliary system, for example an optically pumped defect center [21] or quantum dot [17], which is used to engineer mechanical gain or loss. We assume a general coupling of the form where a is the bosonic operator of the mechanical mode and A is an operator of the auxiliary system. In the absence of the coupling the auxiliary system relaxes with a rate Γ a λ into a steady state described by the density operator ρ 0 a . This separation of timescales allows us to adiabatically eliminate the dynamics of the auxiliary system and derive an effective equation of motion for the mechanical mode. Following the standard treatment of semi-classical laser theory (see, e.g., section 9.3.2 of reference [45]), the result of such a calculation is a Fokker-Planck equation Here P (α) defined via ρ m = d 2 α P (α)|α α| is the Glauber-Sudarshan phase-space representation or P -distribution [29,45] of the mechanical density operator ρ m . In this equation the total gain function and the total diffusion function are defined in terms of the unperturbed expectation value and correlation function of the operator A, where . . . ss = Tr a {. . . ρ 0 a }. The damping rate γ and diffusion rate 2γN th encapsulate the effect of thermal noise from the environment, whereas D q (α) is the intrinsic quantum noise associated with the gain process. Equation (A2) is a standard Fokker-Planck equation and as such it is equivalent to an Itô stochastic differential equation with drift Γ(α) and diffusion D(α) [30, section 4.3.5]. That is, the dynamics of the mechanical system may be modeled using the Itô stochastic differential equation where dW C = 2 −1/2 (dW 1 + idW 2 ) is a complex Wiener increment; dW 1 and dW 2 are independent standard Wiener processes [30, section 3.8.5]. By extending this analysis to two modes with mutual coupling g we then obtain the two coupled equations By formally introducing the white-noise forces F ± (t) = dW C,± /dt, equations (A6) and (A7) reproduce equation (1).

Drift and diffusion for a phonon laser induced by a three-level defect
Reference [21] considers a nitrogen-vacancy (NV) defect centre, which is coupled via strain to the bending motion of a diamond nanobeam. The NV center is modeled as a three-level system with a ground state |g and two near-degenerate electronically excited states |x and |y . If the excited-state splitting ∆E xy matches the vibration frequency of the phonon mode ω m then resonant transitions are induced between the two states. By optically exciting the upper level |x this process excites the phonon mode (i.e., effects gain) whereas optically exciting the lower level |y leads to a cooling process (see Fig. 1 c) in the main text). In the gain configuration the system can be described by the Hamiltonian ( = 1) where λ is the strain coupling constant. For the cooling configuration we obtain the same type of coupling, but with the role of |x and |y reversed. By assuming a radiative decay rate Γ a for both of the two excited states, ω m = ∆E xy (gain configuration) and weak driving Ω Γ a we obtain [21] Γ q (α) = Γ where Γ = 2λ 2 Ω 2 /Γ 3 a and n 0 = Γ 2 a /4λ 2 .
In the following we apply standard dynamical analysis to understand the phases mentioned above and predict the transition points. To do so we first evaluate the possible stationary solutions r ss and z ss of equations (B1) and (B2), which are given by the solutions of gz ss = [Γ(r ss ) − γ]r ss , gr ss = [Γ(z ss ) + γ]z ss .
The stability of these fixed points is then analyzed using the trace-determinant plane of the Jacobian, i.e. the dynamical matrix of the system linearized about said stationary state. The Jacobian for our system is where the prime denotes the derivative. The trace and determinant of J are respectively. Since the eigenvalues of J may be written entirely in terms of τ and δ thus evaluating τ and δ at a particular stationary state fully characterizes its stability and local dynamical structure. For example, if the real part of both λ + and λ − is negative then the stationary state considered is stable, and a non-vanishing imaginary part effects a curl in the phase portrait (the r-z-plane). We summarize this characterization in figure 7 a).

Phases
Phase I: The state r ss = z ss = 0 is an obvious stationary state of equations (B1) and (B2). Substituting this into equations (B6) and (B7) yields (to first order in γ) Consulting figure 7 a), we see that this corresponds to a stable spiral for Γ < g, but for Γ > g becomes a saddle node. We thus conclude that phase I is exhibited for The dashed lines denote nullclines (the curves in the r-z-plane specified by setting eitherṙ = 0 orż = 0; the intersections of nullclines correspond to stationary states), and the solid lines with arrows denote example integral curves (actual trajectories in r and z). c) Cartoon of the general bifurcation diagram for r for 1 < ν < 2. Solid blue lines denote stable stationary states, dashed red lines unstable stationary states, and the shaded blue and red regions denote stable and unstable limit cycles respectively. This plot has been made using data produced by the software MatCont [47] and we have chosen ν = 1.5 and γ/g = 0.001. We note the presence of intermediate phases in the transition from phase IIIw to phase III, the green-shaded region. In this region the stable stationary state is such that r < z, which vitiates our assumption |α| > |β| that we used to eliminate the phase; for |α| < |β| the phase φ = π/2 is unstable. We therefore expect some kind of oscillation in this intermediate regime, which is also observed in numerical simulations.
Phase II: In phase II the two occupation numbers are roughly equal, however simply assuming r ss = z ss yields an inconsistency due to finite γ. Let us therefore substitute the ansatz r ss = r We neglect terms of order γ 2 and higher. Note that, this stationary state only exists for Γ/g ≥ 1.
Substituting equation (B11) into equations (B6) and (B7) yields and We see that δ > 0 only for 1 < Γ/g < [ν/(ν − 1)] ν . However, unlike the previous analysis of phase I, in this case τ can be positive whilst δ is positive. Such is Consulting figure 7 a), one sees that as Γ/g is increased the stationary state (B11) goes from a stable spiral to an unstable spiral and then to a saddle node.
We thus conclude that phase II is exhibited for (B15) Phases IIIw and III: In both phase IIIw and III the occupation numbers are not even roughly equal. To analyze this regime it is instructive to consider the nullclines (the curves in the r-z-plane specified by setting eitherṙ = 0 orż = 0; the intersections of nullclines correspond to stationary states): The function f (r) is always positive, and only zero at r = 0 and r → ∞. The maximum is but the limit cycle persists slightly beyond this upper bound. Since it is only possible that [ν/(ν − 1)] ν > [(ν + 2ν 2 + √ 2ν + 3ν 2 )/(2ν 2 − 1)] ν if 1 ≤ ν < 2, this phase cannot be observed for ν ≥ 2.
The limit cycle does not admit a simple analytic form, however by assuming that it is small and 3. c) Fitting distribution in the lasing regime, for Γ/g = 8.6. Note that the y-axis of the last plot has been scaled by 1/3. centered on stationary state (B11), which is the only unstable spiral for this regime, one may approximate its frequency ω osc by the imaginary part of the eigenvalue of the Jacobian evaluated at equation (B11). One finds ω osc ≈ 2 g 3 (Γ − g)/Γ; this result has been numerically verified for the case ν = 1 (not shown). On the other hand, it is also clear that phase III corresponds to the case of four intersections. One may show that there are four intersections only if Γ/g > [ν/(ν −1)] ν , and it appears that (at least for ν ≥ 2) of the two extra stationary states one must be an unstable spiral and the other a stable spiral. We thus conclude that phase III is exhibited for (B18) Note that this phase cannot be observed for ν = 1. For ν = 2 one may easily check that, given Γ/g > [ν/(ν − 1)] ν = 4, the stable stationary amplitudes are conditions.
From the numerical data we obtain the marginal steady-state distributions P ss (α) and P ss (β) for the two modes. Since the system is invariant under a combined rotation of α and β in phase space, also the marginal distributions are radially symmetric. Figure 8 shows examples of the resulting radial distribution P ss (r = |α|) for the gain mode, before, close to, and after the transition.
To obtain a simple characterization of the system's steady state, we fit the radial distribution P ss (r = |α|) by a distribution of the form P r 0 ,σ (r) = N × r × e − (r−r 0 ) 2 where N is a normalization constant. It corresponds to a thermal state for r 0 = 0 and σ = √ N th and approaches the distribution of a coherent state with random phase in the limit σ r 0 . The optimal values for r 0 and σ are obtained by minimization of the squared difference where H i represents the hight of the i'th bar of the histogram.
where the average |α| 2 is calculated directly from the data set, to represent the range of fluctuations.