${\mathscr{P}}{\mathscr{T}}$-symmetry breaking in the steady state of microscopic gain–loss systems

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 ${ \mathcal P }{ \mathcal T }$ (parity- and time-reversal) symmetry. Their observation triggered considerable interest in discrete and continuous systems with ${ \mathcal P }{ \mathcal T }$-symmetry along with alternative non-Hermitian formulations of quantum theory [2]. Such fundamental considerations remain speculative, but there exist many classical systems in which ${ \mathcal P }{ \mathcal T }$-symmetric dynamics can be obtained with appropriately engineered loss and gain. This was first pointed out in the context of photonic waveguides [3–6], lattices, and resonators [7–10]. Other examples include cold atoms [11–13] and optomechanical devices [14–16]. Of particular interest in such systems is the ${\text{}}{breaking}$ of ${ \mathcal P }{ \mathcal T }$-symmetry, i.e., when by tuning a parameter the energy spectrum becomes complex and the eigenvectors no longer exhibit the underlying ${ \mathcal P }{ \mathcal T }$-symmetry. This phenomenon was first experimentally observed in optical waveguides [5, 6], and is currently the subject of intense experimental and theoretical research.

In this work we go beyond this dynamical picture and address an interesting and still open question: what are the steady states of actual physical systems with ${ \mathcal P }{ \mathcal T }$-symmetry? This question becomes especially important for atomic or microscopic solid-state realizations of gain–loss systems. Here nonlinear saturation effects as well as the presence of thermal and quantum noise have a crucial influence on the system's dynamics and the long-time behavior can no longer be inferred from an eigenvalue analysis only. By focusing on the experimentally relevant example of coupled mechanical resonators with optically induced gain and loss (see figure 1) we show that ${ \mathcal P }{ \mathcal T }$-symmetry breaking in the steady state exhibits various unexpected features and in general occurs via additional intermediate phases with retained or 'weakly' broken ${ \mathcal P }{ \mathcal T }$-symmetry. Most 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 ${ \mathcal P }{ \mathcal T }$-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 ${ \mathcal P }{ \mathcal T }$-symmetric systems operated at low amplitudes and close to the quantum regime.

Zoom In Zoom Out Reset image size

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 ${\omega }_{{\rm{m}}}$ and they are mutually coupled, e.g., mechanically via the support, with strength g. In addition, optical or electrical cooling [17–21] and pumping [21–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 ${\omega }_{{\rm{m}}}$, the semiclassical dynamics of the system is then described by the Itô stochastic differential equation [29, 30]

$$\begin{eqnarray}&&\left(\begin{array}{c}\dot{\alpha }\\ \dot{\beta }\end{array}\right)=\left(\begin{array}{cc}{{\rm{\Gamma }}}_{+}(\alpha ) & -{\rm{i}}g\\ -{\rm{i}}g & {{\rm{\Gamma }}}_{-}(\beta )\end{array}\right)\left(\begin{array}{c}\alpha \\ \beta \end{array}\right)+\left(\begin{array}{c}{F}_{+}(t)\\ {F}_{-}(t)\end{array}\right),\end{eqnarray} \tag{ 1 }$$

where α and β are the dimensionless amplitudes of the pumped and cooled mode respectively (more details on the derivation of equation (1) are given in appendix A).

The optically induced gain and loss rates considered here are of the form

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\pm }(\alpha )=\pm \displaystyle \frac{{\rm{\Gamma }}}{{(1+| \alpha {| }^{2}/{n}_{0})}^{\nu }}-\gamma ,\end{eqnarray} \tag{ 2 }$$

where Γ is the maximal rate and $\sqrt{{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 $\nu =2$ [17]. Instead, for conventional laser amplification with inverted two-level systems we would obtain $\nu =1$ [29]. Finally, γ is the bare mechanical damping rate. Since we are interested in the ${ \mathcal P }{ \mathcal T }$-symmetric limit (defined below), we assume $\gamma /{\rm{\Gamma }}\to 0$. However, in all our calculations we retain a finite $\gamma \gt 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}_{\pm }(t)$ represent two independent white-noise processes with $\langle {F}_{\pm }^{* }(t){F}_{\pm }(t^{\prime} )\rangle ={D}_{\pm }\delta (t-t^{\prime} )$. For resonators coupled to a reservoir of temperature T the diffusion rates are ${D}_{+}(\alpha )={D}_{{\rm{q}}}(\alpha )+2\gamma {N}_{\mathrm{th}}$ and ${D}_{-}=2\gamma {N}_{\mathrm{th}}$, where ${N}_{\mathrm{th}}={({{\rm{e}}}^{{\hslash }{\omega }_{{\rm{m}}}/{k}_{B}T}-1)}^{-1}$. The contribution ${D}_{{\rm{q}}}(\alpha \to 0)=2{\rm{\Gamma }}$ for the gain mode represents the intrinsic quantum noise associated with any amplification process and suggests that noise is a fundamental property of ${ \mathcal P }{ \mathcal T }$-symmetric systems [31–33].

By ignoring the effect of noise and assuming $| \alpha {| }^{2},| \beta {| }^{2}\ll {n}_{0}$, equation (1) can be written as ${\partial }_{t}\psi =-{\rm{i}}H\psi$, with a state vector $\psi ={(\alpha ,\beta )}^{{\rm{T}}}$ and a non-Hermitian 'Hamiltonian'

$$\begin{eqnarray}H=\left(\begin{array}{cc}{\rm{i}}{\rm{\Gamma }} & g\\ g & -{\rm{i}}{\rm{\Gamma }}\end{array}\right).\end{eqnarray} \tag{ 3 }$$

This Hamiltonian is invariant under a combined parity (${ \mathcal P }$) and time-reversal (${ \mathcal T }$) symmetry for all values of ${\rm{\Gamma }}/g$, where ${ \mathcal P }\ :{(\alpha ,\beta )}^{{\rm{T}}}\leftrightarrow {(\beta ,\alpha )}^{{\rm{T}}}$ simply corresponds to an exchange of the two modes and ${ \mathcal T }\ :{\rm{i}}\to -{\rm{i}}$. The eigenvalues and (unnormalized) eigenvectors of H are [34] ${\lambda }_{\pm }=\pm \sqrt{{g}^{2}-{{\rm{\Gamma }}}^{2}}$, ${\psi }_{+}={({{\rm{e}}}^{{\rm{i}}\tfrac{\theta }{2}},{{\rm{e}}}^{-{\rm{i}}\tfrac{\theta }{2}})}^{{\rm{T}}}$ and ${\psi }_{-}={({{\rm{i}}{\rm{e}}}^{-{\rm{i}}\tfrac{\theta }{2}},-{{\rm{ie}}}^{{\rm{i}}\tfrac{\theta }{2}})}^{{\rm{T}}}$, where $\sin (\theta )={\rm{\Gamma }}/g$. We see that despite being non-Hermitian, for ${\rm{\Gamma }}\leqslant g$ both eigenvalues ${\lambda }_{\pm }$ are real, indicating purely coherent oscillations between the two modes. In this regime of unbroken ${ \mathcal P }{ \mathcal T }$-symmetry the eigenvectors are eigenstates of the symmetry operator, i.e., ${ \mathcal P }{ \mathcal T }{\psi }_{\pm }={\psi }_{\pm }$. For ${\rm{\Gamma }}\gt g$, both eigenvalues become imaginary and correspond to a gain and a loss eigenmode. In this regime θ is complex and ${\psi }_{\pm }$ no longer possess the same symmetry as H. One thus speaks of a ${ \mathcal P }{ \mathcal T }$-symmetry-breaking transition occurring at the (exceptional) point ${\rm{\Gamma }}=g$ [2].

This conventionally studied ${ \mathcal P }{ \mathcal T }$-symmetry-breaking effect captures the change in the system's transient dynamics, as 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–38]. Secondly, due to non-decaying oscillations in the ${ \mathcal P }{ \mathcal T }$-symmetric phase, any source of noise would heat up the system indefinitely, and within a linearized model the steady state would again be ill-defined.

Let us now return to the full model in equation (1) and evaluate the stationary state ${\psi }_{\mathrm{ss}}={({\alpha }_{\mathrm{ss}},{\beta }_{\mathrm{ss}})}^{{\rm{T}}}$ in the absence of noise. Figures 2(a) and(b) show the mode occupation numbers $| {\alpha }_{\mathrm{ss}}{| }^{2}$ and $| {\beta }_{\mathrm{ss}}{| }^{2}$ as a function of ${\rm{\Gamma }}/g$, and for the two relevant cases $\nu =2$ and $\nu =1$. Firstly, we observe in both plots the expected transition at ${\rm{\Gamma }}/g{| }_{{\bf{I}}\to {\bf{II}}}=1$. Below this point (phase ${\bf{I}}$) the system dynamics is oscillatory, as expected from the discussion of the ${ \mathcal P }{ \mathcal T }$-symmetric phase above. Note, however, that in contrast to the discussion in section 3 we have now included a finite $\gamma \gt 0$, which leads to a slow overall decay towards a well-defined steady state with $| {\alpha }_{\mathrm{ss}}{| }^{2}=| {\beta }_{\mathrm{ss}}{| }^{2}=0$ (see also figure 2(c)).

Zoom In Zoom Out Reset image size

Above this transition point (phase II) the linearized system dynamics becomes unstable and both resonators reach a finite steady-state occupation number $| {\alpha }_{\mathrm{ss}}{| }^{2}/{n}_{0}=| {\beta }_{\mathrm{ss}}{| }^{2}/{n}_{0}={({\rm{\Gamma }}/g)}^{1/\nu }-1+{ \mathcal O }(\gamma )$, determined by the saturation of ${{\rm{\Gamma }}}_{\pm }(\alpha )$. However, this steady state is still an eigenstate of the symmetry operator, ${ \mathcal P }{ \mathcal T }{\psi }_{\mathrm{ss}}\propto {\psi }_{\mathrm{ss}}$, and contrary to our naïve expectation the system remains ${ \mathcal P }{ \mathcal T }$-symmetric beyond the conventional transition point. The existence of a ${ \mathcal P }{ \mathcal T }$-symmetric steady-state with non-vanishing amplitude can be traced back to the fact that in the present model the ${ \mathcal P }{ \mathcal T }$-symmetry is retained even in the nonlinear regime. This means that for equal amplitudes, $\alpha =\beta$, the gain and the loss rate are still exactly the same, i.e., ${{\rm{\Gamma }}}_{+}(\alpha )=-{{\rm{\Gamma }}}_{-}(\alpha )$, $\forall \alpha$ (see equation (2)). As shown in appendix B, this property of the nonlinear system implies that there exists a symmetric state ${\psi }_{\mathrm{ss}}{({\alpha }_{\mathrm{ss}},{\beta }_{\mathrm{ss}})}^{{\rm{T}}}$, with $| {\alpha }_{\mathrm{ss}}| =| {\beta }_{\mathrm{ss}}|$, which satisfies ${\dot{\psi }}_{\mathrm{ss}}=0$ for all values of ${\rm{\Gamma }}/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 ${ \mathcal P }{ \mathcal T }$-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 $\nu =2$, we find a second transition at ${\rm{\Gamma }}/g=4$, beyond which (phase ${\bf{III}}$)

$$\begin{eqnarray}&&\left.\begin{array}{c}| {\alpha }_{\mathrm{ss}}{| }^{2}\\ | {\beta }_{\mathrm{ss}}{| }^{2}\end{array}\right\}={n}_{0}\times \left(\displaystyle \frac{{\rm{\Gamma }}\pm \sqrt{{\rm{\Gamma }}({\rm{\Gamma }}-4g)}}{2g}-1\right)\end{eqnarray} \tag{ 4 }$$

and the steady state breaks parity, i.e., $| {\alpha }_{\mathrm{ss}}| \ne | {\beta }_{\mathrm{ss}}|$. The ${ \mathcal P }{ \mathcal T }$-symmetry-breaking point can be obtained from a linear stability analysis for phase II and for the present model we find

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\rm{\Gamma }}}{g}\right|}_{{\bf{II}}\to {\bf{III}}}={\left(\displaystyle \frac{\nu }{\nu -1}\right)}^{\nu }.\end{eqnarray} \tag{ 5 }$$

This result would imply the absence of symmetry breaking for $\nu =1$. Instead, for this case we observe a Hopf bifurcation around ${\rm{\Gamma }}/g\approx 5.2$, beyond which the system approaches a limit cycle and undergoes periodic oscillations with the characteristic frequency ${\omega }_{\mathrm{osc}}\approx 2\sqrt{{g}^{3}({\rm{\Gamma }}-g)}/{\rm{\Gamma }}$. The limit cycle is formed symmetrically around a central point ${\psi }_{{\rm{c}}}$ with $| {\alpha }_{{\rm{c}}}| =| {\beta }_{{\rm{c}}}|$, meaning that ${ \mathcal P }{ \mathcal T }$-symmetry 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 (${\bf{IIIw}}$) the ${ \mathcal P }{ \mathcal T }$-symmetry is only 'weakly' broken. In contrast to full symmetry breaking, this transition is related to the finite asymmetry induced by $\gamma \gt 0$, and would be hidden in the idealized models where $\gamma =0$. The value for the transition point is

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\rm{\Gamma }}}{g}\right|}_{{\bf{II}}\to {\bf{IIIw}}}={\left(\displaystyle \frac{\nu +2{\nu }^{2}+\sqrt{2\nu +3{\nu }^{2}}}{2{\nu }^{2}-1}\right)}^{\nu },\end{eqnarray} \tag{ 6 }$$

independent of the precise value of $\gamma \ll {\rm{\Gamma }}$. The general expressions in equations (5) and (6) show that for all intermediate cases, $1\lt \nu \lt 2$, we obtain $1\lt {({\rm{\Gamma }}/g)}_{{\bf{II}}\to {\bf{IIIw}}}\lt {({\rm{\Gamma }}/g)}_{{\bf{II}}\to {\bf{III}}}\lt \infty$ and ${ \mathcal P }{ \mathcal T }$-symmetry breaking occurs in three steps ${\bf{I}}\to {\bf{II}}\to {\bf{IIIw}}\to {\bf{III}}$. More details on the derivation of the above results and the stationary phases for general ν are given in appendix 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 ${ \mathcal P }{ \mathcal T }$-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.

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}_{\pm }=D=2\gamma {N}_{\mathrm{th}}$. For ${\rm{\Gamma }}=0$ thermal noise induces additional amplitude fluctuations of about ${({\rm{\Delta }}\alpha )}^{2}\approx D/(2\gamma )={N}_{\mathrm{th}}$, and we expect that as long as ${N}_{\mathrm{th}}\lt {n}_{0}$ the characteristic features shown in figure 2, which scale with the saturation number n₀, 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}_{\mathrm{th}}\gg {n}_{0}$.

Figure 3(a) shows the results of a numerical simulation of the stochastic equation (1), from which we obtain the steady-state distribution ${P}_{\mathrm{ss}}(\alpha ,\beta )$ for $\nu =2$ and ${N}_{\mathrm{th}}/{n}_{0}=10$. In the following we write $\alpha ={r{\rm{e}}}^{{\rm{i}}{\theta }_{\alpha }}$ and $\beta ={z{\rm{e}}}^{{\rm{i}}{\theta }_{\beta }}$ and make use of the fact that the system dynamics is invariant under a global phase rotation. The exact marginal distributions ${P}_{\mathrm{ss}}(\alpha )$ and ${P}_{\mathrm{ss}}(\beta )$ are then fitted by approximate distributions of the form

$$\begin{eqnarray}&&P(\alpha )\sim {r{\rm{e}}}^{-\tfrac{{(r-{r}_{0})}^{2}}{{\rm{\Delta }}{\alpha }^{2}}},\qquad P(\beta )\sim {z{\rm{e}}}^{-\tfrac{{(z-{z}_{0})}^{2}}{{\rm{\Delta }}{\beta }^{2}}},\end{eqnarray} \tag{ 7 }$$

which allows us to extract a radial shift r₀ and z₀ and the range of fluctuations ${({\rm{\Delta }}\alpha )}^{2}$ and ${({\rm{\Delta }}\beta )}^{2}$ for both modes (see appendix C). From these values plotted in figure 3(a) we see that the thermal noise now completely washes out the features associated with the ${ \mathcal P }{ \mathcal T }$-symmetric phases I and II, and for a large range of Γ the system reaches a steady state (phase ${\bf{T}}$), which is to a good approximation thermal, i.e., ${r}_{0}={z}_{0}=0$ and ${({\rm{\Delta }}\alpha )}^{2}\simeq {({\rm{\Delta }}\beta )}^{2}\simeq {N}_{\mathrm{th}}$. Only after a critical value of ${\rm{\Gamma }}/g{| }_{{\bf{T}}\to {\bf{III}}}\approx 7.5$ are the fluctuations suddenly strongly suppressed. In this regime the system relaxes into an asymmetric coherent state with ${r}_{0}\gt {z}_{0}$ approximately given by the amplitudes $| {\alpha }_{\mathrm{ss}}|$ and $| {\beta }_{\mathrm{ss}}|$ given in equation (4) and ${({\rm{\Delta }}\alpha )}^{2},{({\rm{\Delta }}\beta )}^{2}\sim \gamma {N}_{\mathrm{th}}{\rm{\Gamma }}/{g}^{2},\gamma {N}_{\mathrm{th}}/{\rm{\Gamma }}\ll 1$.

Zoom In Zoom Out Reset image size

Before we proceed, let us connect this transition to the phenomenon of ${ \mathcal P }{ \mathcal T }$-symmetry breaking—hitherto defined only for individual states. To do so we introduce the ${ \mathcal P }{ \mathcal T }$-symmetry parameter

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{\langle {(| \alpha | -| \beta | )}^{2}{\rangle }_{\mathrm{ss}}}{\langle | \alpha {| }^{2}{\rangle }_{\mathrm{ss}}+\langle | \beta {| }^{2}{\rangle }_{\mathrm{ss}}}\leqslant 1,\end{eqnarray} \tag{ 8 }$$

which vanishes for a random set of states ${\psi }_{i}={({\alpha }_{i},{\beta }_{i})}^{{\rm{T}}}$ if and only if each state satisfies ${ \mathcal P }{ \mathcal T }{\psi }_{i}={{\rm{e}}}^{{\rm{i}}{\theta }_{i}}{\psi }_{i}$, with some real phase ${\theta }_{i}$. Figure 3(b) shows that indeed Δ changes at the transition point from ${\rm{\Delta }}\simeq {{\rm{\Delta }}}_{\mathrm{th}}=0.215$ for a thermal state to ${\rm{\Delta }}\to 1$ in the symmetry-broken phase. Note that also in the low-noise limit we obtain ${\rm{\Delta }}={{\rm{\Delta }}}_{\mathrm{th}}\gt 0$ in phase ${\bf{I}}$ and therefore only phase ${\bf{II}}$, where ${\rm{\Delta }}\simeq 0$, has a strictly ${ \mathcal P }{ \mathcal T }$-symmetric steady state.

One of the most striking features visible in figures 3(a) and (c) is that in sharp contrast to a conventional lasing transition, the emerging coherent-state amplitudes after the ${ \mathcal P }{ \mathcal T }$-symmetry breaking point are even smaller than the original level of thermal noise. This surprising effect can be understood as follows. Although at each instant in time the amplitudes $\alpha (t)$ and $\beta (t)$, and therefore the gain and loss rates ${{\rm{\Gamma }}}_{+}(\alpha )$ and ${{\rm{\Gamma }}}_{-}(\beta )$, can be quite different, the average dissipation rate $\bar{{\rm{\Gamma }}}=\langle {{\rm{\Gamma }}}_{-}-{{\rm{\Gamma }}}_{+}{\rangle }_{\mathrm{ss}}$ when evaluated in the thermal phase vanishes, $\bar{{\rm{\Gamma }}}\sim { \mathcal O }(\gamma )\sim 0$. What remains (on average) is the weak coupling to the high-temperature environment. In contrast, in the symmetry-broken phase we have $\langle | \alpha | {\rangle }_{\mathrm{ss}}\gg \langle | \beta | {\rangle }_{\mathrm{ss}}$. Therefore, there is a strong imbalance between loss and gain on average, i.e., $\bar{{\rm{\Gamma }}}\sim { \mathcal O }({\rm{\Gamma }})\gt 0$, and the resulting net cooling effect suppresses fluctuations. Thus, this transition in the average dissipation rate of a stationary system can be seen as the counterpart of the transition from real to imaginary eigenvalues in the conventional definition of ${ \mathcal P }{ \mathcal T }$-symmetry breaking. What we are still missing, however, is a simple criterion, which tells us why the system favors one or the other steady state.

To clarify the mechanism behind the symmetry-breaking transition we focus on the symmetry-broken regime ${\rm{\Gamma }}/g\gg 1$, where we can assume that the amplitude of the gain mode $| \alpha (t)| \approx | {\alpha }_{\mathrm{ss}}|$ and the relative phase $\phi \simeq \pi /2$ are approximately constant. We then obtain the equation of motion for the amplitude of the loss mode⁴ , $z=| \beta |$

$$\begin{eqnarray}&&{\partial }_{t}z=-{\partial }_{z}U(z)+\sqrt{\gamma {N}_{\mathrm{th}}}{\eta }_{z}(t),\end{eqnarray} \tag{ 9 }$$

where $\langle {\eta }_{z}(t){\eta }_{z}(t^{\prime} )\rangle =\delta (t-t^{\prime} )$ and (for $\gamma \to 0$)

$$\begin{eqnarray}&&U(z)=-\displaystyle \frac{{n}_{0}{\rm{\Gamma }}}{2(1+{z}^{2}/{n}_{0})}-g| {\alpha }_{\mathrm{ss}}| z\sin (\phi )-\displaystyle \frac{\gamma {N}_{\mathrm{th}}}{2}\mathrm{log}(z).\end{eqnarray} \tag{ 10 }$$

The function U(z) is an effective potential for z, which is sketched in figure 4(a). This potential has a local minimum at ${z}_{0}=| {\beta }_{\mathrm{ss}}|$ (corresponding to the steady state given in equation (4) for ${N}_{\mathrm{th}}\to 0$), which is separated by a finite barrier ${\rm{\Delta }}U$ from the unstable region $z\gt {z}_{\max }$. In the presence of noise, a system initially located at z ≈ z₀ can escape over this barrier via thermally activated processes with a characteristic rate ${R}_{\mathrm{esc}}\simeq {R}_{0}{{\rm{e}}}^{-\tfrac{2{\rm{\Delta }}U}{\gamma {N}_{\mathrm{th}}}}$, where ${R}_{0}=\sqrt{-U^{\prime\prime} ({z}_{\min })U^{\prime\prime} ({z}_{\max })/(4{\pi }^{2})}$ [30]. This rate increases as Γ is reduced and once ${R}_{\mathrm{esc}}$ exceeds the bare damping γ, any configuration with fixed α and β is rapidly destabilized and the transition to a quasithermal state with strongly fluctuating amplitudes occurs. In figure 4(b) we compare the transition point ${\rm{\Gamma }}/g{| }_{{\bf{T}}\to {\bf{III}}}$ obtained from the condition $\gamma ={R}_{\mathrm{esc}}$ with the numerically evaluated values for various ${N}_{\mathrm{th}}/{n}_{0}\gg 1$. The plot shows that ${ \mathcal P }{ \mathcal T }$-symmetry breaking in the large-noise regime is very well described by this thermal activation model.

Zoom In Zoom Out Reset image size

As pointed out above, the implementation of any gain mechanism is necessarily accompanied by a minimum amount of quantum noise, ${D}_{{\rm{q}}}$, which ensures, e.g., the preservation of Heisenberg's uncertainty relations in a quantum mechanical amplification process. For ${ \mathcal P }{ \mathcal T }$-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}_{\mathrm{th}}=0$ and consider the quantum diffusion term

$$\begin{eqnarray}&&{D}_{{\rm{q}}}(\alpha )=\displaystyle \frac{2{\rm{\Gamma }}}{{(1+| \alpha {| }^{2}/{n}_{0})}^{3}},\end{eqnarray} \tag{ 11 }$$

only. The assumed form of ${D}_{{\rm{q}}}(\alpha )$ applies for the three-level gain mechanism shown in figure 1, i.e., for the case $\nu =2$ (see appendix 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₀, z₀ and the range of fluctuations ${\rm{\Delta }}\alpha$ and ${\rm{\Delta }}\beta$ 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 ${ \mathcal P }{ \mathcal T }$-symmetric phase and a ${ \mathcal P }{ \mathcal T }$-symmetry-broken phase with strongly reduced fluctuations. Thus, we conclude that the ${ \mathcal P }{ \mathcal T }$-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}\gg 1$ while for ${n}_{0}\sim 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.

Zoom In Zoom Out Reset image size

Finally, to show that ${ \mathcal P }{ \mathcal T }$-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 [41, 42], as depicted in figure 6(a). In this case the amplitudes ${\alpha }_{n}$ and ${\beta }_{n}$ for each unit cell obey

$$\begin{eqnarray}&&\left(\begin{array}{c}{\dot{\alpha }}_{n}\\ {\dot{\beta }}_{n}\end{array}\right)=\left(\begin{array}{cc}{{\rm{\Gamma }}}_{+}({\alpha }_{n}) & -{\rm{i}}g\\ -{\rm{i}}g & {{\rm{\Gamma }}}_{-}({\beta }_{n})\end{array}\right)\left(\begin{array}{c}{\alpha }_{n}\\ {\beta }_{n}\end{array}\right)-{{\rm{i}}g}^{\prime }\left(\begin{array}{c}{\beta }_{n-1}\\ {\alpha }_{n+1}\end{array}\right)+\left(\begin{array}{c}{F}_{n,+}(t)\\ {F}_{n,-}(t)\end{array}\right),\end{eqnarray} \tag{ 12 }$$

where $g^{\prime}$ is the coupling between the unit cells and ${F}_{n,\pm }(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 $\nu =2$. The observed features can be understood from the plane wave ansatz [41, 43] ${\alpha }_{n}={A}_{k}{{\rm{e}}}^{{\rm{i}}kn}$, ${\beta }_{n}={B}_{k}{{\rm{e}}}^{{\rm{i}}kn}$, where $k=4\pi 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\mapsto {g}_{k}=| g+g^{\prime} {{\rm{e}}}^{{\rm{i}}k}|$. By considering now each of the modes separately, we see that all the transition points identified above occur first for the mode with the largest value of ratio ${\rm{\Gamma }}/| {g}_{k}|$, or equivalently, the smallest value of $| {g}_{k}|$. Since in our model g and $g^{\prime}$ are assumed to be real and positive, this minimal coupling is alway achieved for the $k=\pi$ mode, which therefore determines the symmetry-breaking properties of the array [42, 44]. For an even number of unit cells and the special case $g=g^{\prime}$, i.e., g_k = 0, the gain and the loss modes completely decouple, and the ${ \mathcal P }{ \mathcal T }$-symmetry breaking transitions, both in the linear and nonlinear regime, already occur at ${\rm{\Gamma }}{| }_{{\bf{I}}\to {\bf{II}}}\simeq {\rm{\Gamma }}{| }_{{\bf{II}}\to {\bf{III}}}\simeq {\rm{\Gamma }}{| }_{{\bf{T}}\to {\bf{III}}}\simeq 0$. For linear ${ \mathcal P }{ \mathcal T }$-symmetry chains this instability at ${\rm{\Gamma }}=0$ has already been pointed out in previous works [9, 42, 44]. 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^{\prime}$. We see that the single-mode ansatz captures well the relevant physics both in the low- and high-noise regime. Note, however, that 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.

Zoom In Zoom Out Reset image size

We have analyzed the breaking of ${ \mathcal P }{ \mathcal T }$-symmetry in the steady state of coupled mechanical systems with loss and gain. Our study shows how saturation effects and the influence of thermal or quantum noise lead to surprising new effects, which are not expected from the conventional eigenvalue analysis of idealized ${ \mathcal P }{ \mathcal T }$-symmetric systems. These findings constitute a significant and important first step towards a general understanding of the behavior of realistic systems with ${ \mathcal P }{ \mathcal T }$-symmetry, especially for systems operated close to the quantum regime or for systems with strong nonlinearities.

The present analysis was mainly motivated by recent developments in the field of opto- and nanomechanical systems, where optical cooling and lasing schemes for mechanical oscillators are currently explored experimentally. In particular, for cooling schemes based on single emitters [17, 18, 21–23], the competition between saturation and thermal noise, i.e., ${n}_{0}\sim {N}_{\mathrm{th}}\sim {10}^{2}\mbox{--}{10}^{4}$, is typically met, which makes these systems a natural platform for studying the symmetry breaking both in the low-noise and high-noise regime. However, our findings are quite general and there are many other atomic and solid state systems, where artificial loss and gain can be engineered. For example, based on the already existing excellent level of control, the use of coupled vibrational modes of trapped ion systems [25, 45] would allow a large degree of tunability and such systems would also be well suited to explore the quantum noise limit [46].

We thank T E Lee, G J Milburn, and H E Türeci for stimulating discussions. This work was supported by the European FP7/ITC Project SIQS (600645), by the People Program (Marie Curie Actions) of the European Union under REA grant agreement number PIOF-GA-2011-303228 (project NOLACOME), project OPSOQI (316607) of the WWTF, the Austrian Science Fund (FWF) through SFB FOQUS F40, SFB NextLite F49-10, project GePartWave No. I 1142-N27, DK CoQuS W 1210 and the START grant Y 591-N16. KGM is also supported by the European Union Seventh Framework Programme (FP7-REGPOT-2012-2013-1) under grant agreement 316165.

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 [17] or quantum dot [18], which is used to engineer mechanical gain or loss. We assume a general coupling of the form

$$\begin{eqnarray}&&{H}_{\mathrm{int}}=\lambda ({{Aa}}^{\dagger }+{A}^{\dagger }a),\end{eqnarray} \tag{ A.1 }$$

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 ${{\rm{\Gamma }}}_{{\rm{a}}}\gg \lambda$ into a steady state described by the density operator ${\rho }_{{\rm{a}}}^{0}$. 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 [47]), the result of such a calculation is a Fokker–Planck equation of the form

$$\begin{eqnarray}&&\dot{P}(\alpha ,t)=\left[-\left(\displaystyle \frac{\partial }{\partial \alpha }\alpha +\displaystyle \frac{\partial }{\partial {\alpha }^{* }}{\alpha }^{* }\right){\rm{\Gamma }}(\alpha )+\displaystyle \frac{{\partial }^{2}}{\partial {\alpha }^{* }\partial \alpha }D(\alpha )\right]P(\alpha ,t).\end{eqnarray} \tag{ A.2 }$$

Here $P(\alpha )$ defined via ${\rho }_{{\rm{m}}}=\int {{\rm{d}}}^{2}\alpha \,P(\alpha )| \alpha \rangle \langle \alpha |$ is the Glauber–Sudarshan phase–space representation or P-distribution [29, 47] of the mechanical density operator ${\rho }_{{\rm{m}}}$. In this equation the total gain function

$$\begin{eqnarray}&&{\rm{\Gamma }}(\alpha )={{\rm{\Gamma }}}_{{\rm{q}}}(\alpha )-\gamma ,\qquad {{\rm{\Gamma }}}_{{\rm{q}}}(\alpha )=-\displaystyle \frac{{\rm{i}}\lambda \langle A{\rangle }_{\mathrm{ss}}}{\alpha }\end{eqnarray} \tag{ A.3 }$$

and the total diffusion function

$$\begin{eqnarray}&&D(\alpha )=2\gamma {N}_{\mathrm{th}}+{D}_{{\rm{q}}}(\alpha ),\qquad {D}_{{\rm{q}}}(\alpha )=2{\lambda }^{2}{\mathfrak{R}}{\int }_{0}^{\infty }{\rm{d}}\tau (\langle {A}^{\dagger }(\tau )A{\rangle }_{\mathrm{ss}}-| \langle A{\rangle }_{\mathrm{ss}}{| }^{2}),\end{eqnarray} \tag{ A.4 }$$

are defined in terms of the unperturbed expectation value and correlation function of the operator A, where $\langle ...{\rangle }_{\mathrm{ss}}={\mathrm{Tr}}_{{\rm{a}}}\{...{\rho }_{{\rm{a}}}^{0}\}$. The damping rate γ and diffusion rate $2\gamma {N}_{\mathrm{th}}$ encapsulate the effect of thermal noise from the environment, whereas ${D}_{{\rm{q}}}(\alpha )$ is the intrinsic quantum noise associated with the gain process.

Equation (A.2) is a standard Fokker–Planck equation and as such it is equivalent to an Itô stochastic differential equation with drift ${\rm{\Gamma }}(\alpha )$ and diffusion $D(\alpha )$ [30, section 4.3.5]. That is, the dynamics of the mechanical system may be modeled using the Itô stochastic differential equation thus

$$\begin{eqnarray}&&{\rm{d}}\alpha ={\rm{\Gamma }}(\alpha )\alpha {\rm{d}}t+\sqrt{D(\alpha )}{{\rm{d}}W}_{{\mathbb{C}}},\end{eqnarray} \tag{ A.5 }$$

where ${{\rm{d}}W}_{{\mathbb{C}}}={2}^{-1/2}({{\rm{d}}W}_{1}+{{\rm{i}}{\rm{d}}W}_{2})$ is a complex Wiener increment; dW₁ and dW₂ 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

$$\begin{eqnarray}&&{\rm{d}}\alpha ={{\rm{\Gamma }}}_{+}(\alpha )\alpha {\rm{d}}t-{\rm{i}}g\beta {\rm{d}}t+\sqrt{{D}_{+}(\alpha )}{{\rm{d}}W}_{{\mathbb{C}},+},\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}&&{\rm{d}}\beta ={{\rm{\Gamma }}}_{-}(\beta )\beta {\rm{d}}t-{\rm{i}}g\alpha {\rm{d}}t+\sqrt{{D}_{-}(\beta )}{{\rm{d}}W}_{{\mathbb{C}},-}.\end{eqnarray} \tag{ A.7 }$$

By formally introducing the white-noise forces ${F}_{\pm }(t)={{\rm{d}}W}_{{\mathbb{C}}},\pm /{\rm{d}}t$, equations (A.6) and (A.7) reproduce equation (1).

A.1. Drift and diffusion for a phonon laser induced by a three-level defect

Reference [17] considers a nitrogen-vacancy (NV) defect center, 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\rangle$ and two near-degenerate electronically excited states $| x\rangle$ and $| y\rangle$. If the excited-state splitting ${\rm{\Delta }}{E}_{{xy}}$ matches the vibration frequency of the phonon mode ${\omega }_{{\rm{m}}}$ then resonant transitions are induced between the two states. By optically exciting the upper level $| x\rangle$ this process excites the phonon mode (i.e., effects gain) whereas optically exciting the lower level $| y\rangle$ leads to a cooling process (see figure 1(c) in the main text). In the gain configuration the system can be described by the Hamiltonian $({\hslash }=1)$

$$\begin{eqnarray}&&H={\omega }_{m}{a}^{\dagger }a-{\rm{\Delta }}{E}_{{xy}}| y\rangle \langle y| +\displaystyle \frac{{\rm{\Omega }}}{2}(| g\rangle \langle x| +| x\rangle \langle g| )+\lambda (| y\rangle \langle x| {a}^{\dagger }+a| x\rangle \langle y| ),\end{eqnarray} \tag{ A.8 }$$

where λ is the strain coupling constant. For the cooling configuration we obtain the same type of coupling, but with the role of $| x\rangle$ and $| y\rangle$ reversed. By assuming a radiative decay rate ${{\rm{\Gamma }}}_{{\rm{a}}}$ for both of the two excited states, ${\omega }_{{\rm{m}}}={\rm{\Delta }}{E}_{{xy}}$ (gain configuration) and weak driving ${\rm{\Omega }}\ll {{\rm{\Gamma }}}_{{\rm{a}}}$ we obtain [17]

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{q}}}(\alpha )=\displaystyle \frac{{\rm{\Gamma }}}{{(1+| \alpha {| }^{2}/{n}_{0})}^{2}},\qquad {D}_{{\rm{q}}}(\alpha )=\displaystyle \frac{2{\rm{\Gamma }}}{{(1+| \alpha {| }^{2}/{n}_{0})}^{3}},\end{eqnarray} \tag{ A.9 }$$

where ${\rm{\Gamma }}=2{\lambda }^{2}{{\rm{\Omega }}}^{2}/{{\rm{\Gamma }}}_{{\rm{a}}}^{3}$ and ${n}_{0}={{\rm{\Gamma }}}_{{\rm{a}}}^{2}/4{\lambda }^{2}$.

In this section we outline the stability analysis used to obtain the steady-state phases of two coupled resonators in the absence of noise, i.e., with ${D}_{\pm }=0$. By moving to polar coordinates $\alpha ={r{\rm{e}}}^{{\rm{i}}{\theta }_{\alpha }}$ and $\beta ={z{\rm{e}}}^{{\rm{i}}{\theta }_{\beta }}$, equation (1) can be rewritten as

$$\begin{eqnarray}&&\dot{r}=-[\gamma -{\rm{\Gamma }}(r)]r-g\sin (\phi )z,\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray}&&\dot{z}=-[\gamma +{\rm{\Gamma }}(z)]z+g\sin (\phi )r,\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&\dot{\phi }=g\left(\displaystyle \frac{r}{z}-\displaystyle \frac{z}{r}\right)\cos (\phi ),\end{eqnarray} \tag{ B.3 }$$

where $\phi ={\theta }_{\alpha }-{\theta }_{\beta }$. Note that the system is invariant under a combined rotation of α and β and therefore the evolution of the total phase ${\theta }_{\alpha }+{\theta }_{\beta }$ can be neglected. For the last equation we see that there are two fixed points for the phase, ${\phi }_{\mathrm{ss}}=\pm \pi /2$. Due to finite γ, the stationary occupation number of the gain mode $| \alpha {| }^{2}$ is always slightly larger than that of the loss mode $| \beta {| }^{2}$, wherefore ${\phi }_{\mathrm{ss}}=\pi /2$ is the stable. We therefore set $\phi =\pi /2$ and henceforth study the two-dimensional dynamical system with variables r and z.

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}_{\mathrm{ss}}$ and ${z}_{\mathrm{ss}}$ of equations (B.1) and (B.2), which are given by the solutions of

$$\begin{eqnarray}&&{{gz}}_{\mathrm{ss}}=[{\rm{\Gamma }}({r}_{\mathrm{ss}})-\gamma ]{r}_{\mathrm{ss}},\qquad {{gr}}_{\mathrm{ss}}=[{\rm{\Gamma }}({z}_{\mathrm{ss}})+\gamma ]{z}_{\mathrm{ss}}.\end{eqnarray} \tag{ B.4 }$$

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

$$\begin{eqnarray}{\bf{J}}=\left(\begin{array}{cc}{\rm{\Gamma }}^{\prime} ({r}_{\mathrm{ss}}){r}_{\mathrm{ss}}+{\rm{\Gamma }}({r}_{\mathrm{ss}})-\gamma & -g\\ g & -{\rm{\Gamma }}^{\prime} ({z}_{\mathrm{ss}}){z}_{\mathrm{ss}}-{\rm{\Gamma }}({z}_{\mathrm{ss}})-\gamma \end{array}\right),\end{eqnarray} \tag{ B.5 }$$

where the prime denotes the derivative. The trace and determinant of ${\bf{J}}$ are

$$\begin{eqnarray}&&\tau =\mathrm{Tr}{\bf{J}}=-2\gamma +{\rm{\Gamma }}^{\prime} ({r}_{\mathrm{ss}}){r}_{\mathrm{ss}}-{\rm{\Gamma }}^{\prime} ({z}_{\mathrm{ss}}){z}_{\mathrm{ss}}+{\rm{\Gamma }}({r}_{\mathrm{ss}})-{\rm{\Gamma }}({z}_{\mathrm{ss}}),\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&\delta =\det {\bf{J}}={g}^{2}-[{\rm{\Gamma }}^{\prime} ({r}_{\mathrm{ss}}){r}_{\mathrm{ss}}+{\rm{\Gamma }}({r}_{\mathrm{ss}})-\gamma ][{\rm{\Gamma }}^{\prime} ({z}_{\mathrm{ss}}){z}_{\mathrm{ss}}+{\rm{\Gamma }}({z}_{\mathrm{ss}})+\gamma ]\end{eqnarray} \tag{ B.7 }$$

respectively. Since the eigenvalues of ${\bf{J}}$ may be written entirely in terms of τ and δ thus

$$\begin{eqnarray}&&{\lambda }_{\pm }=\tfrac{1}{2}(\tau \pm \sqrt{{\tau }^{2}-4\delta }),\end{eqnarray} \tag{ B.8 }$$

evaluating τ and δ at a particular stationary state fully characterizes its stability and local dynamical structure. For example, if the real part of both ${\lambda }_{+}$ and ${\lambda }_{-}$ 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 B1(a).

Zoom In Zoom Out Reset image size

B.1. Phases

Phase I:

The state ${r}_{\mathrm{ss}}={z}_{\mathrm{ss}}=0$ is an obvious stationary state of equations (B.1) and (B.2). Substituting this into equations (B.6) and (B.7) yields (to first order in γ)

$$\begin{eqnarray}&&\tau =-2\gamma ,\qquad \delta ={g}^{2}-{{\rm{\Gamma }}}^{2}.\end{eqnarray} \tag{ B.9 }$$

Consulting figure B1(a), we see that this corresponds to a stable spiral for ${\rm{\Gamma }}\lt g$, but for ${\rm{\Gamma }}\gt g$ becomes a saddle node. We thus conclude that phase I is exhibited for

$$\begin{eqnarray}&&0\lt \displaystyle \frac{{\rm{\Gamma }}}{g}\lt 1.\end{eqnarray} \tag{ B.10 }$$

Phase II:

In phase II the two occupation numbers are roughly equal, however simply assuming ${r}_{\mathrm{ss}}={z}_{\mathrm{ss}}$ yields an inconsistency due to finite γ. Let us therefore substitute the ansatz ${r}_{\mathrm{ss}}={r}_{\mathrm{ss}}^{(0)}+\gamma {r}_{\mathrm{ss}}^{(1)}+\cdots$ and ${z}_{\mathrm{ss}}={r}_{\mathrm{ss}}^{(0)}+\gamma {r}_{\mathrm{ss}}^{(1)}+\cdots$ and equate alike-orders in γ. The result is

$$\begin{eqnarray}{r}_{\mathrm{ss}}^{(0)} & = & \sqrt{{n}_{0}}\times \sqrt{{({\rm{\Gamma }}/g)}^{1/\nu }-1},\\ {r}_{\mathrm{ss}}^{(1)} & = & -{z}_{\mathrm{ss}}^{(1)}={r}_{\mathrm{ss}}^{(0)}{[{\rm{\Gamma }}^{\prime} ({r}_{\mathrm{ss}}^{(0)}){r}_{\mathrm{ss}}^{(0)}+{\rm{\Gamma }}({r}_{\mathrm{ss}}^{(0)})+g]}^{-1}.\end{eqnarray} \tag{ B.11 }$$

We neglect terms of order ${\gamma }^{2}$ and higher. Note that, this stationary state only exists for ${\rm{\Gamma }}/g\geqslant 1$. Substituting equation (B.11) into equations (B.6) and (B.7) yields

$$\begin{eqnarray}&&\tau =-2\gamma \left(1-\displaystyle \frac{{\rm{\Gamma }}^{\prime\prime} ({r}_{\mathrm{ss}}^{(0)}){r}_{\mathrm{ss}}^{(0)2}+2{\rm{\Gamma }}^{\prime} ({r}_{\mathrm{ss}}^{(0)}){r}_{\mathrm{ss}}^{(0)}}{{\rm{\Gamma }}^{\prime} ({r}_{\mathrm{ss}}^{(0)}){r}_{\mathrm{ss}}^{(0)}+{\rm{\Gamma }}({r}_{\mathrm{ss}}^{(0)})+g}\right)\end{eqnarray} \tag{ B.12 }$$

$$\begin{eqnarray}&&=-2\gamma \left(1-\displaystyle \frac{2\nu [{({\rm{\Gamma }}/g)}^{1/\nu }-1][(2\nu -1){({\rm{\Gamma }}/g)}^{1/\nu }-2(\nu +1)]}{2\nu {({\rm{\Gamma }}/g)}^{1/\nu }+2(1-\nu ){({\rm{\Gamma }}/g)}^{2/\nu }}\right)\end{eqnarray} \tag{ B.13 }$$

and

$$\begin{eqnarray}&&\delta ={g}^{2}-{[{\rm{\Gamma }}^{\prime} ({r}_{\mathrm{ss}}^{(0)}){r}_{\mathrm{ss}}^{(0)}+{\rm{\Gamma }}({r}_{\mathrm{ss}}^{(0)})]}^{2}={g}^{2}-{g}^{2}{[2\nu -1-2\nu {({\rm{\Gamma }}/g)}^{-1/\nu }]}^{2}.\end{eqnarray} \tag{ B.14 }$$

We see that $\delta \gt 0$ only for $1\lt {\rm{\Gamma }}/g\lt {[\nu /(\nu -1)]}^{\nu }$. However, unlike the previous analysis of phase I, in this case τ can be positive while δ is positive. Such is the case if ${[\nu /(\nu -1)]}^{\nu }\gt {[(\nu +2{\nu }^{2}+\sqrt{2\nu +3{\nu }^{2}})/(2{\nu }^{2}-1)]}^{\nu }$. Consulting figure B1(a), one sees that as ${\rm{\Gamma }}/g$ is increased the stationary state (B.11) goes from a stable spiral to an unstable spiral and then to a saddle node. We thus conclude that phase II is exhibited for

$$\begin{eqnarray}&&1\lt \displaystyle \frac{{\rm{\Gamma }}}{g}\lt \min \left\{{\left(\displaystyle \frac{\nu +2{\nu }^{2}+\sqrt{2\nu +3{\nu }^{2}}}{2{\nu }^{2}-1}\right)}^{\nu },{\left(\displaystyle \frac{\nu }{\nu -1}\right)}^{\nu }\right\}.\end{eqnarray} \tag{ B.15 }$$

Phases ${\bf{IIIw}}$ and ${\bf{III}}$:

In both phase ${\bf{IIIw}}$ and ${\bf{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 $\dot{r}=0$ or $\dot{z}=0;$ the intersections of nullclines correspond to stationary states):

$$\begin{eqnarray}&&z=f(r)-\displaystyle \frac{\gamma }{g}r,\qquad r=f(z)+\displaystyle \frac{\gamma }{g}z,\qquad f(r)=\displaystyle \frac{{\rm{\Gamma }}}{g}\displaystyle \frac{r}{{(1+{r}^{2}/{n}_{0})}^{\nu }}.\end{eqnarray} \tag{ B.16 }$$

The function f(r) is always positive, and only zero at r = 0 and $r\to \infty$. The maximum is located at $r/\sqrt{{n}_{0}}={(2\nu -1)}^{-1/2}$ where $f{(r)={({\rm{\Gamma }}/g)(2\nu )}^{-1}(2\nu -1)}^{\nu +1/2}$. The additional terms $-(\gamma /g)r$ and $(\gamma /g)z$ approximately rotate the r- and z-axes for $\gamma \ll g$. The possible intersections of the two nullclines as a function of ${\rm{\Gamma }}/g$ are shown in figure B1(c): one intersection (I in the figure), two intersections (II and ${\bf{IIIw}}$ in the figure), or four intersections (III in the figure). It is clear that phase ${\bf{IIIw}}$ corresponds to the case of two intersections, and from the previous section we know that if ${[(\nu +2{\nu }^{2}-\sqrt{2\nu +3{\nu }^{2}})/(2{\nu }^{2}-1)]}^{\nu }\lt {\rm{\Gamma }}/g\lt {[\nu /(\nu -1)]}^{\nu }$ then neither of these two intersections corresponds to a stable stationary state. Therefore, by the Poincaré–Bendixon theorem [48, section 7.3], a limit cycle must exist. We thus conclude that phase ${\bf{IIIw}}$ is exhibited for

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\nu +2{\nu }^{2}+\sqrt{2\nu +3{\nu }^{2}}}{2{\nu }^{2}-1}\right)}^{\nu }\lt \displaystyle \frac{{\rm{\Gamma }}}{g}\lt {\left(\displaystyle \frac{\nu }{\nu -1}\right)}^{\nu },\end{eqnarray} \tag{ B.17 }$$

but the limit cycle persists slightly beyond this upper bound. Since it is only possible that ${[\nu /(\nu -1)]}^{\nu }\gt {[(\nu +2{\nu }^{2}+\sqrt{2\nu +3{\nu }^{2}})/(2{\nu }^{2}-1)]}^{\nu }$ if $1\leqslant \nu \lt 2$, this phase cannot be observed for $\nu \geqslant 2$. The limit cycle does not admit a simple analytic form, however by assuming that it is small and centered on stationary state (B.11), which is the only unstable spiral for this regime, one may approximate its frequency ${\omega }_{\mathrm{osc}}$ by the imaginary part of the eigenvalue of the Jacobian evaluated at equation (B.11). One finds ${\omega }_{\mathrm{osc}}\approx 2\sqrt{{g}^{3}({\rm{\Gamma }}-g)}/{\rm{\Gamma }};$ this result has been numerically verified for the case $\nu =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 ${\rm{\Gamma }}/g\gt {[\nu /(\nu -1)]}^{\nu }$, and it appears that (at least for $\nu \geqslant 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

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Gamma }}}{g}\gt {\left(\displaystyle \frac{\nu }{\nu -1}\right)}^{\nu }.\end{eqnarray} \tag{ B.18 }$$

Note that this phase cannot be observed for $\nu =1$. For $\nu =2$ one may easily check that, given ${\rm{\Gamma }}/g\gt {[\nu /(\nu -1)]}^{\nu }=4$, the stable stationary amplitudes are

$$\begin{eqnarray}&&{r}_{\mathrm{ss}}/\sqrt{{n}_{0}}={\{[{\rm{\Gamma }}+\sqrt{{\rm{\Gamma }}({\rm{\Gamma }}-4g)}]/(2g)-1\}}^{1/2},\end{eqnarray} \tag{ B.19 }$$

$$\begin{eqnarray}&&{z}_{\mathrm{ss}}/\sqrt{{n}_{0}}={\{[{\rm{\Gamma }}-\sqrt{{\rm{\Gamma }}({\rm{\Gamma }}-4g)}]/(2g)-1\}}^{1/2}.\end{eqnarray} \tag{ B.20 }$$

For the numerical results presented in figures 3 and 4 in the main text we have simulated the stochastic equations (A.6) and (A.7) (and the corresponding extended set of equations for the array) using the Euler–Maruyama method (see section 15 of [30]). Since we are interested in the steady state of the system, we collect data (complex numbers $\alpha (t)$ and $\beta (t)$) after a time of $5\times {\tau }_{0}$, where ${\tau }_{0}=1/\gamma$ is approximately the time scale in which the steady state is approached. Over a period of $45\times {\tau }_{0}$ we select 4000 random data points. We repeat the procedure 80 times with random initial conditions.

From the numerical data we obtain the marginal steady-state distributions ${P}_{\mathrm{ss}}(\alpha )$ and ${P}_{\mathrm{ss}}(\beta )$ 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 C1 shows examples of the resulting radial distribution ${P}_{\mathrm{ss}}(r=| \alpha | )$ for the gain mode, before, close to, and after the transition.

Zoom In Zoom Out Reset image size

To obtain a simple characterization of the system's steady state, we fit the radial distribution ${P}_{\mathrm{ss}}(r=| \alpha | )$ by a distribution of the form

$$\begin{eqnarray}&&{P}_{{r}_{0},\sigma }(r)={ \mathcal N }\times r\times {{\rm{e}}}^{-\tfrac{{(r-{r}_{0})}^{2}}{{\sigma }^{2}}},\end{eqnarray} \tag{ C.1 }$$

where ${ \mathcal N }$ is a normalization constant. It corresponds to a thermal state for ${r}_{0}=0$ and $\sigma =\sqrt{{N}_{\mathrm{th}}}$ and approaches the distribution of a coherent state with random phase in the limit $\sigma \ll {r}_{0}$. The optimal values for r₀ and σ are obtained by minimization of the squared difference

$$\begin{eqnarray}&&{\chi }^{2}({r}_{0},\sigma )=\displaystyle \sum _{i=1}^{N}{({H}_{i}-{P}_{{r}_{0},\sigma }({r}_{i}))}^{2},\end{eqnarray} \tag{ C.2 }$$

where H_i represents the hight of the ith bar of the histogram.

While this fitting procedure yields accurate values for the radial displacement of the distribution maximum, r₀, we find that the corresponding values for the width σ do not very well capture the broad thermal background in the vicinity of the transition (see figure C1(b)). Therefore, we use instead the quantity

$$\begin{eqnarray}&&{({\rm{\Delta }}\alpha )}^{2}=\langle | \alpha {| }^{2}\rangle -{r}_{0}^{2},\end{eqnarray} \tag{ C.3 }$$

where the average $\langle | \alpha {| }^{2}\rangle$ is calculated directly from the data set, to represent the range of fluctuations.