Topological dynamics and excitations in lasers and condensates with saturable gain or loss.

We classify symmetry-protected and symmetry-breaking dynamical solutions for nonlinear saturable bosonic systems that display a non-hermitian charge-conjugation symmetry, as realized in a series of recent groundbreaking experiments with lasers and exciton polaritons. In particular, we show that these systems support stable symmetry-protected modes that mirror the concept of zero-modes in topological quantum systems, as well as symmetry-protected power-oscillations with no counterpart in the linear case. In analogy to topological phases in linear systems, the number and nature of symmetry-protected solutions can change. The spectral degeneracies signalling phase transitions in linear counterparts extend to bifurcations in the nonlinear context. As bifurcations relate to qualitative changes in the linear stability against changes of the initial conditions, the symmetry-protected solutions and phase transitions can also be characterized by topological excitations, which set them apart from symmetry-breaking solutions. The stipulated symmetry appears naturally when one introduces nonlinear gain or loss into spectrally symmetric bosonic systems, as we illustrate for one-dimensional topological laser arrays with saturable gain and two-dimensional flat-band polariton condensates with density-dependent loss.

A wide range of topological quantum effects manifest themselves in symmetries of an excitation spectrum. This relation has been explored most extensively for fermionic systems, for which the single-particle Hamiltonian can obey a chiral symmetry or a charge-conjugation symmetry [1][2][3]. Allowing also for the possibility of time-reversal invariance one can identify ten universality classes featuring a variety of topological quantum numbers. These quantum numbers determine the formation of zero modes and unidirectional transport channels at edges, surfaces, interfaces and defects of various dimensions [4,5]. In the fermionic case these symmetries are fundamental as they are owed to the structure of Fock or Nambu space [6], and therefore also extend to the interacting case [7] . A prime example are superconducting systems, for which the charge-conjugation symmetry of the Bogoliubov-De Gennes equation dicates that excitations ψ exp(−iωt) at a positive frequency ω are paired with excitations Xψ * exp(iωt) at the corresponding negative frequency −ω, where X is a suitable unitary transformation. The system can then also feature robust Majorana zero modes ψ = Xψ * at ω = 0 [8][9][10], an effect which translates to the fermion parity anomaly in the full many-body theory [11][12][13].
Some of these symmetries also play a natural role beyond the fermionic context. For instance, time-reversal symmetry is intimately related to optical reciprocity [14], a classical concept whose modification leads to topological effects in photonic structures [15] and analogous optical [16][17][18], acoustic [19,20] and mechanical [21] systems. Topological effects can also be engineered into bosonic quantum systems, such as cold atomic gases [22] and exciton polaritons [23][24][25]. A variant of the chargeconjugation symmetry indeed applies to the excitation spectrum of interacting quantum systems with a complex classical-wave limit, such as in the mean-field analysis of weakly interacting bosons. The excitations are again characterized by a Bogoliubov theory [26][27][28][29], and recent work has explored a variety of mechanisms to engineer topological features into this description [30][31][32][33][34][35][36]. Furthermore, exploiting the concept of a pseudospin, the charge-conjugation symmetry can also be induced into the single-particle theory of bosonic systems, where one again can admit linear gain or loss [37][38][39][40][41]. The spectrum of the effective Hamiltonian then displays symmetric pairs of complex frequencies (Ω n , −Ω * n ), which can be interpreted as resonances. Topological features persist because this spectral constraint can enforce a number of purely imaginary self-symmetric resonances Ω n = −Ω * n , which provide an analogue to broadened fermionic zero modes [42][43][44][45].
Here, we consider bosonic quantum systems for which a notion of charge-conjugation symmetry applies directly to the time-dependent nonlinear wave equation. Nonlinearities are an often unavoidable feature of systems with gain and loss, and the question whether topological states retain their properties in such settings is both fundamentally and practically important. Extending the class of systems with a pseudospin to feature nonlinear gain or loss, we are naturally led to consider nonlinear complex-wave dynamics where the time-dependent solutions appear in pairs Ψ(t) and XΨ * (t). In contrast to a general time-reversal symmetry, which induces a partner solution XΨ * (−t), the solutions Ψ(t) and XΨ * (t) both exhibit the same arrow of time. This timepreserving symmetry allows for self-symmetric statesincluding stationary states Ψ(t) = Ψ 0 = XΨ * 0 , which we interpret as zero modes, and time-dependent oscillating states with Ψ(t + T ) = Ψ(t) = XΨ * (t)-as well as a twisted variant of periodically oscillating states with Ψ(t + T /2) = XΨ * (t). As we will see, the stability of these solutions and their competition with the symmetry-breaking pairs is governed by topological excitations that allow us to define different phases. In particular, the Bogoliubov theory of stationary zero modes features an additional charge-conjugation symmetry that renders the system topologically nontrivial and results in a constrained excitation spectrum. A similar constraint occurs for self-symmetric oscillating solutions, while their twisted variant features topologically protected excitations that arise from the propagation over a half-period. We describe these concepts in general terms and illus- Im ω Re ω Re ω trate our findings for two model systems, representing one-dimensional topological laser arrays with saturable gain and two-dimensional flat-band polariton condensates with density-dependent loss (see Figs. 1

and 2).
General model and classification of states.-Consider a bosonic system whose classical limit is described by a complex-wave equation The operator H provides an effective single-particle description, while the functional F describes the nonlinear effects which may encompass gain and loss. We assume these effects to be local in time and often will drop the time argument. In keeping with the quantum origin of the wave function, the nonlinear effects should preserve the global U(1) gauge freedom, so that any solution can be multiplied by a fixed phase factor exp(iα). This can be guaranteed by assuming F Ψ Ψ = F Ψ * Ψ * , where subscripts denote functional derivatives. In addition, we here stipulate that H displays a charge-conjugation symmetry, XHX = −H * where X is a unitary operator with X 2 = 1. To extend this notion to the nonlinear case, we similarly demand that XF [Ψ, Ψ * ]X = −(F [XΨ * , XΨ]) * , with the same operator X.
We will describe practical ways to realize such systems soon below but first discuss the consequences of the following general observation. For any solution Ψ(t) of Eq. (1), we obtain a partner solutioñ Given our assumptions this correspondence even applies when H and F contain an explicit time dependence [46]; however, we focus on the autonomous case. Let us first reflect on the possible classes of solutions admitted by Eq. (1). Because of the underlying charge-conjugation symmetry of H, the linear system possesses pairs of solutions Ψ(t) = Ψ n exp(−iΩ n t) and Ψ(t) = XΨ * n exp(iΩ * n t), where the latter expression corresponds to a frequencyΩ n = −Ω * n . This includes the possibility of topologically protected zero modes with purely imaginary frequency Ω n =Ω n = i Im Ω n and time dependence Ψ(t) =Ψ(t) = Ψ n exp(Im Ω n t). If for any of these states Im Ω n is positive the linear system is unstable, but the nonlinear effects can stabilize the system. This results in stationary states Ψ(t) = Ψ n exp(−iΩ n t) with real Ω n , which are self-consistent solutions of the condition [47] The solutions either (a) occur in pairs Ψ n ,Ψ n where Ω n andΩ n = −Ω n now are both real, or (b) are time-independent zero modes Ψ 0 =Ψ 0 whose frequency Ω 0 = 0 now vanishes [48]. Alternatively, the system may tend to time-dependent solutions, including periodic, aperiodic and chaotic solutions. These will either (c) still occur in pairs Ψ(t),Ψ(t) that bear no further relation besides Eq. (2), or may be constrained in two possible ways. (d) The time-dependent solution may be self-symmetric, Ψ(t) =Ψ(t); given this condition at some point of time, it will remain true for all times -0.6 0 -5 5 Im ω Re ω Re ω [49]. (e) Two partner solutions may be related by a finite time-offset, Ψ(t) =Ψ(t + T /2). It then follows that the solutions must be periodic, Ψ(t + T ) = Ψ(t), which amounts to the twisted states mentioned in the introduction. The nature of a state can be assessed, e.g., by comparing the correlation functions C = | Ψ(0)|Ψ(t) |, C = | Ψ(0)|Ψ(t) |, which coincide or alternate for selfsymmetric or twisted solutions. Figures 1 and 2 illustrate these different types of solutions for two model systems, which are both formed on bipartite lattices (N A sites A and N B sites B) that give rise to a pseudospin degree of freedom. The dynamics is governed by the equations where we grouped the amplitudes on both sublattices into vectors A, B. The real nearest-neighbour couplings t kl form a linear Hamiltonian H that supports at least ν = |N A − N B | zero modes, irrespective of whether the system is homogenous or not [50].
In the examples this describes a finite segment of a Su-Schrieffer-Heeger (SSH) chain with alternating couplings t k,k+1 = 1, 0.7 (Fig. 1), where a topological midgap state arises from a central coupling defect [39,51,52], or a finite segment of a disordered Lieb lattice (Fig. 2), which exhibits multiple zero modes associated with a flat band [53,54]. Both of these models can be realized on a large variety of platforms, including bosonic systems [40,[55][56][57][58][59][60][61][62][63][64][65][66][67]. For the SSH chain, we consider nonlinearities that represent saturable gain, as often adopted in laser models, while for the Lieb lattice we consider density-dependent loss, as considered, e.g., in studies of polaritonic condensates. For any solution with amplitudes A(t) and B(t), Eq. (4) exhibits a partner solution with amplitudes A * (t) and −B * (t), so that X = 1 1 A ⊕ (−1 1 B ) acts as a Pauli-z matrix in pseudospin space. Admitting different amounts of gain g A , g B and loss γ A , γ B on the two sublattices allows us to study the competition between topological and nontopological states, which leads to the different examples shown in the figures.
Note that the form of X in these models implies that self-symmetric states display a rigid phase-shift of i between the two sublattices. Given thatΨ(t) = Ψ(t), the amplitudes A(t) are real while the amplitudes B(t) are imaginary, which then persists for all times. On the other hand, symmetry-breaking stationary states with a finite frequency must have equal weight on both sublattices [39]. The different types of states can therefore also be discriminated by their distinct spatial features.
Topological features of the excitation spectrum.-To complete the classification of these solutions we consider their stability, which can be addressed by applying Bogoliubov theory [68]. This amounts to identifying the eigenmodes of linear perturbations, and results in a spectrum of 2(N A + N B ) excitation frequencies ω m . For a stable system all of these excitations obey Im ω m ≤ 0.
For non-self-symmetric pairs of stationary modes Ψ n , Ψ n , the time-preserving symmetry of Eq. (1) implies that their excitation spectra are identical, but does not impose any further constraints on these spectra. For selfsymmetric stationary zero modes Ψ 0 =Ψ 0 with Ω 0 = 0, on the other hand, we can classify the perturbations into joint eigenstates of X with eigenvalue σ = ±1. For our models of gain or loss, the perturbations fulfill the reduced eigenvalue equations where f is diagonal with entries f A (|A k | 2 ) or f B (|B k | 2 ) (depending on the sublattice), and the prime denotes the derivative with respect to the argument. The excitation spectrum is thus composed of two parts, eigenvalues ω + that display an enhanced nonlinearity and eigenvalues ω − that coincide with those of the nonlinear wave operator H + f , with Ψ = Ψ 0 fixed to the zero mode. Panel (b) in Figs. 1 and 2 shows that the eigenvalues ω − determine the stability of the stationary zero modes while the eigenvalues ω + lie deep in the complex plane. Therefore, if we restrict the perturbations to preserve the symmetry, the stability of such modes is greatly enhanced. Note that both reduced spectra in these examples have an odd number of excitations. Thus, each reduced spectrum has an odd number of eigenvalues pinned to the imaginary axis. Setting u − = Ψ 0 we recover that one of the eigenvalues ω −,0 = Ω 0 = 0 vanishes, in accordance with the U(1) symmetry of the wave equation.
A similar analysis can be carried out for timedependent solutions Ψ(t), whose stability is then characterized by a corresponding quasienergy spectrum obtained from a propagator U (t) with eigenvalues λ m = exp(−iω m t). This again includes the U(1) Goldstone mode with λ 0 = 1, but also an additional Goldstone mode λ T = 1 that describes the time-translation freedom of any solution Ψ(t). For general pairs of states Ψ(t),Ψ(t), the time-preserving symmetry again dictates that their quasienergy excitation spectra are identical. For self-symmetric states Ψ(t) =Ψ(t), we can once more separate perturbations that preserve or break the symmetry. This leads to time-dependent variants of Eq. (6), where f and f inherit their time dependence from Ψ(t). The stability spectrum thus again contains a component ω − inherited from the nonlinear wave operator H + f , Im ω Re ω twisted-mode side which includes the vanishing eigenvalue ω −,0 = 0 (with u −,0 = Ψ(t)) dictated by the U(1) gauge freedom. In addition, there is now also an eigenvalue ω +,0 ≡ ω T = 0 arising from the time-translation freedom of any solution Ψ(t), which is associated with u +,0 = dΨ(t)/dt.
For the final case of a twisted periodic state with Ψ(t + T /2) =Ψ(t), the evolution operator of the perturbations factorizes as U (T ) = U 2 with a half-step propagator U = X Σ x U (T /2). The U(1) gauge and timetranslation freedoms can now in principle be satisfied in two ways, corresponding to ω 0,T T /2 = 0 or π. Inspecting the Goldstone mode ψ 0 (t) = (Ψ(t), −Ψ * (t)) T we invariably find ψ 0 (T /2) = −X Σ x ψ 0 (0), so that the twisted solutions always realize the case ω 0 T /2 = π, while the timetranslation freedom remains associated with ω T T /2 = 0. This separation of eigenvalues to symmetry-protected positions is again a robust spectral feature that can only change in topological phase transitions, an example of which is illustrated in Fig. 3.
Conclusions.-In summary, we described how the charge-conjugation symmetry can be generalized to nonlinear complex-wave equations, where it can be interpreted as an antilinear symmetry that preserves the arrow of time. This provides a topological interpretation of a range of modes, including stationary and oscillating self-symmetric solutions as well as twisted oscillating modes, which display a characteristically constrained stability spectrum that governs their competition with symmetry-breaking solutions. These properties allow us to identify topological states in weakly interacting bosonic systems with saturable gain or densitydependent loss, such as realized in topological lasers or polaritonic condensates in structures that exhibit a suitable pseudospin degree of freedom. In such models the self-symmetric states exhibit a rigid phase relation between different wave components, which also equips them with a distinct spatial feature. These findings show that the notion of symmetry-protected topological states persists in nonlinear bosonic systems, and that the nonlinearities lead to an even richer phenomenology than in the original single-particle context. We gratefully acknowledge support by EPSRC via Programme Grant No. EP/N031776/1 and Grant No. EP/P010180/1.
Here we provide a more detailed derivation of the stability excitation spectrum for the different solutions identified in the main text.
For stationary states, their stability can be analyzed by setting (8) and expanding in the small quantities u and v, which we group into a vector ψ = (u, v) T . This leads to the Bogoliubov equation where Σ z = σ z ⊗ 1 1 is a Pauli matrix in the space of u and v while the Bogoliubov Hamiltonian is given by In our models with saturable gain or density-dependent loss [Eq. (5) in the main text], this becomes  The model based on the Lieb lattice also consists of 21 sites, but these are arranged in two dimensions so that 12 are A sites and 9 are B sites. In the hermitian limit, there are now at least three zero modes, even in the case of disorder in the couplings. In this model we study density-dependent losses, and include disorder in the couplings t kl as well as in the gain and loss parameters gA,B and γA,B. This disorder is generated from independent random numbers rn with a box distribution over [0.5, 1.5], so that pn = prn for any model parameter pn with average p.
where the terms containing f have to be read as diagonal matrix with entries f A (|A k | 2 ) or f B (|B k | 2 ) (depending on the sublattice), and the prime denotes the derivative with respect to the argument. As any Bogoliubov equation, Eq. (9) exhibits the charge-conjugation symmetry with the corresponding Pauli matrix Σ x = σ x ⊗ 1 1. This implies that the excitation spectrum is symmetric under the transformation ω m →ω m = −ω * m . Furthermore, due to the U(1) gauge freedom, Eq. (9) always admits the Goldstone mode ψ 0 = (Ψ, −Ψ * ) T with eigenvalue pinned at ω 0 = 0. As the dimension of H is even, an additional odd number of eigenvalues must lie on the imaginary axis. We need to explore the interplay of these spectral features with the analogous symmetries in Ω.
The time-preserving symmetry of the nonlinear wave equation (1) (main text) induces the relation where X = 1 1 2 ⊗ X. For non-self-symmetric pairs of stationary modes Ψ n ,Ψ n , this implies that their excitation spectra are identical. For self-symmetric stationary zero modes Ψ 0 =Ψ 0 with Ω 0 = 0, on the other hand, Eq. (11) turns into an additional charge-conjugation symmetry which imposes a constraint on the excitation spectrum.
In this case, we can classify the perturbations into joint eigenstates of X with eigenvalue σ = ±1. These take the form v σ = σXu σ , and fulfill the reduced eigenvalue equations For our models with saturable gain or densitydependent loss, this simplifies to Eq. (6) in the main text. As described there, the excitation spectrum is thus composed of two parts, eigenvalues ω + that display an enhanced nonlinearity and eigenvalues ω − that coincide with those of the nonlinear wave operator H + f , with Ψ = Ψ 0 fixed to the zero mode. According to Eq. (8), the perturbations [e −iω−t u − − e iω−t Xu * − ] corresponding to ω − break the time-preserving symmetry of the zero mode. Therefore, if we restrict the perturbations to preserve the symmetry, the stability of such modes is greatly enhanced. Setting u − = Ψ 0 we recover that one of the eigenvalues ω − = Ω 0 = 0 vanishes, in accordance with the U(1) symmetry of the wave equation.
A similar analysis can be carried out for timedependent solutions Ψ(t). Their stability against perturbations u(t) + v * (t) is governed by the time-dependent Bogoliubov equation whose solutions ψ(t) = U (t)ψ(0) define the propagator U (t). The eigenvalues λ m = exp(−iω m t) of U (t) can be represented by complex quasienergies that are constrained in similar ways as the excitations in the stationary case . The conventional charge-conjugation symmetry H * (t) = −Σ x HΣ x implies U * = Σ x U Σ x , so that each ω m comes with a partnerω m = −ω * m (thusλ m = λ * m ) or is purely imaginary (whereupon λ m is real). This includes the U(1) Goldstone mode ψ 0 (t) = (Ψ(t), −Ψ * (t)) T , as well as an additional Goldstone mode ψ T (t) = (dΨ/dt, dΨ * /dt) T which describes the time-translation freedom t 0 of any solution Ψ(t + t 0 ).
For general pairs of states Ψ(t),Ψ(t), the propagators are related as U (t) = XŨ * (t)X = X Σ xŨ (t)Σ x X , which again implies that their quasienergy excitation spectra are identical. This also applies to the case of periodic solutions with Ψ(t) = Ψ(t + T ), where the propagator U (T ) represents a Bogoliubov-Floquet operator. Such solutions are unstable if |λ m | < 1, hence Im ω m > 0, apart from the two Goldstone modes who are fixed at λ 0 = λ T = 1. For the more general case of periodic solutions with Ψ(t) = exp(iϕ)Ψ(t + T ), this applies when the Floquet operator is redefined as diag(e −iϕ , e iϕ )U (T ).
For self-symmetric states Ψ(t) =Ψ(t), we have at all times U (t) = X U * (t)X = X Σ x U (t)Σ x X . This allows us to separate perturbations that preserve or break the symmetry, and leads to the time-dependent variant (7) of Eq. (6), both given in the main text. For periodic states Ψ(t) = exp(iϕ)Ψ(t + T ), the self-symmetry constraints the phase to ϕ = 0, π, hence τ = exp(iϕ) = ±1, which needs to be respected in the proper definition of the Bogoliubov-Floquet operator F = τ U (T ) to result in two Goldstone modes pinned at λ 0 = λ T = 1.