Non-Hermitian Topological End-Mode Lasing in Polariton Systems

We predict the existence of non-Hermitian topologically protected end states in a one-dimensional exciton-polariton condensate lattice. The topological transitions are driven by spatially modulated incoherent pump which imposes a four-site unit cell. We show that in particular limits of the model the number of end states is described by a Chern number and a topological invariant based on the eigenvalues of the Wilson loop. By continuity the end states exist beyond these limits because the end modes can only disappear or appear at transitions where the energy gap between the bulk bands and end states closes. We find that such transitions in an infinite open chain arise due to enforced exceptional points which can be predicted directly from the bulk Bloch wave functions. This allows us to establish a new type of bulk-boundary correspondence for non-Hermitian systems and to compute the phase diagram of an infinite open chain analytically. Finally, we demonstrate topological lasing of a single end-mode in a realistic model of a microcavity lattice.

In this Letter we study NH physics in a light-matter hybrid quantum quasiparticle system [36] and demonstrate the feasibility of NH topological end states and end mode lasing. The non-trivial topology is achieved through spatial modulation of external incoherent pump intensity in a one-dimensional lattice of homogeneously coupled microcavity exciton-polaritons. These are natural candidate systems for realizing topological phases [37,38], but in existing proposals [39][40][41][42][43][44][45][46][47][48][49][50][51][52] and realizations [28,[53][54][55][56] topological order arises from the bandtopology in a Hermitian Hamiltonian, so it is not related to the open dissipative NH nature of the system. We show that the presence of excitonic component opens up a possibility to realize a polariton system with competing Hermitian and non-Hermitian effects, so that unique topological phases possessing various number of topological end states can be observed. In particular limits of the model the number of end states is described by a Chern number [23] and a topological invariant based on the eigenvalues of the Wilson loop [57], but by continuity the end states exist beyond these special cases. We find that the end modes can disappear or appear at transitions where the energy gap between the bulk bands and the end states closes without closing of the bulk energy gap. Such transitions can occur in open NH systems due to the NH skin effect [17][18][19]22]. However, we find that skin effect is not present in our model because the Hamiltonian satisfies a specific type of NH time-reversal symmetry [15]. We discover that the transitions occurring in an infinite open chain arise because specific bulk Bloch wave functions are incompatible with the boundary conditions leading to enforced exceptional points in the spectrum. This way we establish a new type of bulkboundary correspondence for NH Hamiltonians, because the changes in the number of end states can be directly predicted from the bulk Bloch wave functions. This powerful new tool allows us to compute the topological phase diagram of an infinite open system analytically.
Moreover, we demonstrate the possibility of lasing of end modes for realistic physical parameters corresponding to a lattice of polariton micropillars. By appropriate choice of pumping intensity pattern, we achieve the situation when only the selected end state is amplified, which results in topological single mode lasing.
Theoretical model. We consider a one-dimensional (1D) lattice of coupled micropillars as depicted in Fig. 1(a). Each micropillar contains a quantum well and is assumed to host a tightly bound exciton-polariton mode [58]. We start with the system of discrete mean-arXiv:1912.06431v1 [cond-mat.mes-hall] 13 Dec 2019 In the special cases where the bulk bands and the end states can be ordered based on the real part of the energies we use a notation (m1, m2, m3), where mi is the number of end states in the ith energy gap. (d),(e) Numerically and analytically computed phase diagrams as a function of g1 and g2 for θ = π/3. The small disagreement between the phase diagrams arises due to the finite size effects. All energies are in units of κ.
field Gross-Pitaevskii equations [28,59] i ψ n = −κ nn ψ m + g c |ψ n | 2 + g R n R n + i where ψ n (t) is the condensate amplitude in the n-th lattice cell, the nn sum runs over nearest neighbors, n R n (t) is the density of exciton reservoir in the n-th cell, P n is the external nonresonant pumping rate, γ c and γ R are the decay rates of the condensate and the reservoir, respectively, g c and g R are the corresponding interaction constants, and R is the rate of scattering from the reservoir to the condensate. We assume that the polariton interactions within the condensate represented by the term g c |ψ n | 2 are negligible in comparison with the reservoircondensate interaction g R n R n , which is a good approximation in most experiments where nonresonant pumping is used.
While the coupling coefficients κ are assumed to be the same for each pair of neighboring micropillars, the external pumping P n is modulated spatially with the four-site periodicity. Therefore, the pumping gives rise to complex on-site potential [ Fig. 1(a)]. This leads to an effective NH tight-binding Hamiltonian [60,61] Here k is the Bloch wave number in the units of inverse intercell distance, 1,2 = g 1,2 e iθ are the complex on-site potentials, and g n and θ satisfy The real part of 1,2 stands for the interaction-induced on-site potential, whereas the imaginary part corresponds to the balance between gain and loss. In the linear approximation of negligibly small |ψ n | 2 , which we assume in most of this work, the reservoir occupation n R n in the above equation is defined by the local pumping rate and takes the form n R n = n R n0 = P n /γ R , and we have chosen the pump rates P n so that the onsite potentials in every second lattice sites are related to each other. The real amplitudes g n can be positive or negative depending on the local pumping intensity [60]. For realistic repulsive interactions between polaritons one has 0 < θ < π/2. When complex on-site potentials are absent ( 1,2 = 0), the lattice is trivial since all the coupling coefficients κ are equal. In typical polariton micropillar lattices κ ∼ 0.1 meV.
Topological phase diagram. We have numerically implemented the Hamiltonian (2) with open boundary conditions and calculated the number of end states as a function of various parameters of the model [60]. The results are shown in Figs on the real part of the energies we also use a notation (m 1 , m 2 , m 3 ), where m i is the number of end states in the energy gap between ith and i + 1th bands.
In the extreme NH limit θ = π/2 the system obeys NH chiral symmetry SĤ k (g 1 , g 2 )S = −Ĥ † k (g 1 , g 2 ) with S = 1⊗σ z . In this case the number of end states with zero real part of the energy is described by a Chern number C [23,60]. For g 2 > 0 and g 1 = |g 2 | we obtain C = −1 which means that there should be one state with zero real part of the energy at each end of the chain in agreement with our numerical calculations [Figs. 1(c) and 2(a)]. Thus, this limiting case of the model belongs to the [1+1] phase with the number of end states in each energy gap given by (0, 2, 0). In this limit the system exhibits topological states solely due to the non-Hermiticity [34].
If θ = 0, g 1 = |g 2 | and g 2 > 0 the system is Hermitian and obeys an inversion symmetry IĤ k (g, g) Therefore, a topological invariant can be defined based on a Wilson loop invariant or quantized Zak phases [57,60]. We find that both approaches yield identical predictions [60]. Namely, we find that the first two energy gaps are trivial whereas the third energy gap is nontrivial, so that we expect one topological mode at each end of the chain in agreement with our numerical calculations [ Fig. 2(b)]. Thus, this limiting case of the model belongs to [1+1] phase with the number of end states in each energy gap given by (0, 0, 2) [ Fig. 1(c)]. In this case the system exhibits topological states due to the Hermitian topological invariants.
If θ = 0, g 1 = |g 2 | and g 2 < 0 the bulk Hamiltonian obeys inversion symmetry but the unit cell is not compatible with the inversion symmetry. This means that in the case of open boundary conditions one needs to include additionally half of a unit cell in the end of the chain in order to make the Hamiltonian inversion symmetric. This obscures the bulk-boundary correspondence but our non-rigorous reasoning based on the Wilson loop invariants indicates that this limit of the model belongs to [2+2] phase with the number of end states in each energy gap given by (1, 2, 1) [60]. These invariants are in agreement with our numerical calculations [ Fig. 2(c)].
Although in a realistic description of exciton-polaritons with a pump-induced potential both the real and the imaginary parts of the potential are non-negligible (0 < θ < π/2), the limits considered above are interesting because the end states can only disappear or appear at transitions where the energy gap between the bulk bands and end states closes. Therefore, we immediately obtain the overall structure of the phases [1+1] and [2+2] shown in Fig. 1(b),(c). There is an energy gap closing between bands if g 1 = g 2 = 0, so this phase transition line separating the [1+1] and [2+2] phases is also easy to understand.
We find that if π/4 < θ < π/2 there exists additional transitions where the end modes disappear or appear without closing of the bulk gap [ Fig. 1(b),(c)]. Such transitions are known to occur due to the NH skin effect [17][18][19]22], but the skin effect is absent in our model because of NH time-reversal symmetryĤ * −k (g 1 , g 2 ) =Ĥ † k (g 1 , g 2 ) [15]. Instead, we notice that the eigenstates of the continuum bands in an infinite open chain can be constructed from standing waves of the bulk Bloch wave functions. This construction usually follows without any problems along the lines suggested in Ref. [22] but we discover that if the first or last component of the bulk Bloch wave function vanishes it is impossible to satisfy the boundary conditions. Therefore, states should disappear from the continuum bands at the values of the model parameters where such incompatibility with the boundary conditions occurs. This can happen in NH systems with the help of exceptional points which transform into new end states. In the following we demonstrate that with the help of this reasoning we can compute the topological phase diagram of an infinite open system analytically.
To find the enforced exceptional points in an infinite open chain we need to calculate when the product Q k of the first and the last components of all eigenvectors vanishes. In the case g 1 = |g 2 | we find by a straightforward calculation that Therefore, Q k vanishes along the contour The vanishing Q k leads to appearance of two exceptional points because of inversion (chiral-inversion) symmetry if g 2 > 0 (g 2 < 0). Therefore, the number of end states changes by 2 at these transition lines. By computing the end state spectrum for phase when g 2 and θ are tuned across the transition line defined by Eq. (4) [ Fig. 1(b),(c) and Fig. 2(d),(e)]. In addition to the prediction of the phase boundaries, we can predict the momentum and the band where the exceptional point leads to appearance or disappearance of the end states. We have checked that the numerical calculations are in agreement with these analytical considerations.
Similarly we can find the zeros of Q k in the g 1 -g 2 plane for fixed θ. The explicit expression for Q k is much more complicated but we arrive to a simple result that Q k vanishes along the four elliptic contours defined by and g 1 = ±κ tan θ sin k 2 , g 2 = ∓κ cos k 2 , k ∈ [−π, π).
In the general case g 1 = |g 2 | the vanishing of Q k leads to a single exceptional point so that the number of end states changes by 1 at these transition lines. Topological lasing. We now consider the possibility of lasing of the NH end states. We expect that single-mode lasing can be realized if all modes decay in time except a selected end mode, which is amplified by the medium. This requires an adjustment of the spectra presented in Fig. 2. As demonstrated in Supplementary Information [60], the following set of coupled equations describes evolution in the weakly nonlinear regime i ψ n = −κ nn ψ m + n ψ n −Γ n (1 + i tan θ) |ψ n | 2 ψ n , (5) where now we include additional loss γ in each node, i. e. A = 1 − iγ, B = − 2 − iγ, C = − 1 − iγ, and This results in the shift of the imaginary part of the excitation spectrum by −γ. For a careful choice of γ, we can stabilize all the eigenstates except one. Physically, it corresponds to homogeneous reduction of pumping across the lattice. The reduction of pumping affects also the real part of the potential, but this can be eliminated by introducing a rotating frame for the condensate amplitudes [60]. The parameter Γ n = P n Rg R /γ 2 R describes the nonlinear interactions with the reservoir, which includes both the nonlinear decay and repulsive interparticle interactions.
In Fig. 3 we show examples of the evolution of the intensity distributions obtained by solving Eq. (5) in a L = 4 × 10 unit-cell configuration with a random initial amplitude. After initial evolution, the polariton density saturates at a steady-state distribution, which approximately corresponds to the shape of the corresponding end mode in the linear spectrum, see Fig. 3(c).
In the symmetric case with g 1 = |g 2 | the end states are always below the bulk states in the spectra of imaginary parts. Hence, to establish single-end mode lasing, one has to apply additional pump to one of the end pillars [28]. This leads to the increase of the imaginary part of the corresponding end state, relative to all other states, and can lead to its single mode amplification. The resulting imaginary part of the spectrum is shown in the inset of 3(b). Here the values of on-site effective potentials g = ±0.14meV can be created with external pumps P + ≈ 7.9meVµm −2 ps −1 , and P − ≈ 2.75meVµm −2 ps −1 .
In addition, we apply extra pump to the rightmost pillar with δP ≈ 3.6meVµm −2 ps −1 .
Conclusions. We predict the emergence of NH topological states in one-dimensional arrays of coupled polariton micropillars as a result of periodic spatial modulation of the incoherent pump rate. We have established the topological phase diagram

SUPPLEMENTARY MATERIAL
Derivation from the polariton model.-The system of coupled micropillars with polariton resonance is described by Eq. (1) of the main text. In the approximation of negligibly small polariton interactions, the first of the two equations can be expressed as The prefactor of last term is equivalent to the on-site effective energy of the tight-binding Hamiltonian (2). One has The last term in the above is a constant real energy shift, which can be removed by introducing a rotating frame for condensate amplitudes, ψ n → ψ n e −i(g R γc/R)t . In result, we obtain the effective complex potentials in the form with tan θ = R 2g R . It is clear that while g 1,2 can be either positive or negative, the parameter θ must fulfill tan θ > 0. In the linear approximation of negligibly small |ψ n | 2 , n R n in the above equation is simply proportional to the local pumping rate and takes the form It should be noted that the Eq. (6) can be obtained from the effective tight-binding Hamiltonian whereâ n ,b n ,ĉ n andd n (â † n ,b † n ,ĉ † n andd † n ) denote annihilation (creation) operators for the sub-lattices A, B, C and D in the nth lattice cell of our structure, 1 ≤ n ≤ L.
Here L is the number of unit cells, 4L is the total number of sites, and i are the complex effective potentials. By going to the momentum space we arrive to the Eq. (2) in the main text. To go beyond the linear approximation, one may employ the adiabatic approximation [61] where (11) in the first order approximation with respect to the condensate density.
Numerical procedure for end mode detection.-In the main text we discuss presence of topologically protected end modes in the case of open boundary conditions [ Fig. 1(b),(d) and Fig. 2 in main text]. For each values of the parameters (g 1 , g 2 , θ), we have first calculated eigenvalues and eigenstates by solving the eigenvalue problem for the Hamiltonian (10), and then classified the eigenstates to bulk and end modes by evaluating the ratio between the integral of the intensity in the most far left (and most far right) unit cell with the total eigenmodes norm. For bulk (end) states this ratio goes to zero (nonzero) in the limit L → ∞, but for finite L we identify the states as end states if the ratio is above a chosen threshold. We have used L = 100 and threshold ratios 0.04 and 0.2 in Fig. 1(b) and Fig. 1(d), respectively.
Chern number.-In the case θ = π/2 the Hamiltonian (2) in the main text satisfies a NH chiral symmetry with S = 1 ⊗ σ z . This means that we can define a bulk topological invariant which describes the number of end modes [23]. Namely, we can obtain an effective two-dimensional Hermitian Hamiltonian as and in the case of open boundary conditions the zeroenergy end states of H eff (k, η) correspond to end states ofĤ k (g 1 , g 2 ) with zero real part of energy and imaginary part equal to η [23]. Because H eff (k, η) is defined in a two-dimensional (k, η) space and its Chern number C is quantized if it is gapped, C = 0 at half-filling implies that the HamiltonianĤ k (g 1 , g 2 ) supports |C| end states with zero real part of the energy [23]. The Chern number C for H eff (k, η) can be obtained using Kubo formula where the Berry curvature Ω k,η is given by Here ψ p k,η are eigenstates of H eff (k, η) (sorted in ascending order of eigenenergy) and n F is the number of occupied bands.
Wilson loop invariant.-If g 1 = |g 2 | the bulk Hamiltonian obeys an inversion symmetry I. We focus here on the case g 2 > 0. In this case IĤ k (g, g)I −1 =Ĥ −k (g, g), where I = σ x ⊗ σ x . Therefore, we can define a topological invariant N (−1) using the ideas of Ref. [57]. The most general definition makes use of Wilson loop operator given by the matrix elements where i, j = 1, . . . , n F label the occupied bands ψ (i) k of H k , and P k is the projector on occupied bands for a given k. The step size δk has to be chosen sufficiently small. For non-orthogonal eigenvectors projector P k can be easily constructed after orthogonalizing the set of occupied eigenvectors that span the same subspace so According to [57] in the Hermitian case and in the presence of inversion symmetry the spectrum of W consists of eigenvalues that are either ±1 or pairs of complex conjugate numbers λ and λ . The topological invariant N (−1) is then given by the number of −1 eigenvalues of W. We have checked that the spectra of W hold the same properties in the NH case and hence we can define the N (−1) invariant in an analogous way. An important simplification coming from [57] is that the invariant can be expressed as where n (−1) (0) and n (−1) (π) are the numbers of occupied states at k = 0 and k = π having inversion eigenvalue −1 (I commutes with H k at these points). We have checked that in our case this formula also gives identical results as the Wilson loop approach. We note that in order to use these invariants in description of non-Hermitian systems one needs to be able to keep track which bands are occupied as a function of the parameters of the model. In the main text we use these invariants only in the Hermitian case where the bands can be ordered based on the energies. This way we can define a topological invariant for each energy gap by assuming that all bands below the gap are occupied when calculating N (−1) . By collecting the invariants calculated for the different gaps we obtain a Wilson loop invariant {N (−1),1 , N (−1),2 , N (−1),3 }.
The approach can be easily generalized to all situations where the bands can be ordered according to the real (or imaginary) parts of the energies. In a more general situation, one needs to be able to follow the spectral flow of the bands as a function of the parameters of the model.
Now, assuming inversion symmetry and δk → 0 one can show that proving the quantization of Zak phase and providing simplification with respect to original formula. The Zak phases of the four bands seem to predict correctly the (non)triviality of a given energy gap in the sense that when the bands can be ordered based on the real part of the energy we have the following correspondence between Wilson loop invariants of the energy gaps and Zak phases of the bands In terms of the Zak phases the energy gap is topologically nontrivial if the product of the Zak phases of the occupied bands is −1.
Point and line gap closing.-If one does not keep track of the spectral flow of the occupied bands in the calculation of the Zak phases there is an important difference between non-Hermitian and Hermitian systems. Namely, in Hermitian systems the Zak phases can only change when two bands touch each other, but in non-Hermitian systems the Zak phases can change in two different ways: (i) bands touch in the complex-energy plane (point gap closing) and exchange their invariants or (ii) bands are interchanged in the complex plane without touching but with closing of a line energy gap (Fig. 4). If one follows the spectral flow of the bands as a function of the parameters of the model and keeps the same bands occupied the Zak phases can only change in a point gap closing. If g 1 = |g 2 | and g 2 > 0 we find that m i = 2N (−1),i (i = 1, 2, 3) because each nontrivial Wilson loop invariant gives rise to topological state at both ends of the chain. In the case of our model we obtain N (−1),1 = N (−1),2 = 0 and N (−1),3 = 1 so that the number of end modes in each energy gap are (0, 0, 2).
However, the situation is slightly more complicated if g 1 = |g 2 | and g 2 < 0. In this case the bulk Hamiltonian 0 Im E ν 2 =+1 ν 3 =-1 phase (0,0,1) phase (0,1,1) phase (0,0,1) phase (0,1,1) obeys inversion symmetry but the unit cell is not compatible with the inversion symmetry. This means that in the case of open boundary conditions one needs to include additionally half of a unit cell in the end of the chain in order to make the Hamiltonian inversion symmetric. This obscures the bulk-boundary correspondence because one of the ends comes from a bulk Hamiltonian where parameters are replaced as g i → −g i . Numerically, we find that for θ = 0, g 1 = |g 2 | and g 2 < 0 the Wilson loop invariants are {1, 1, 0} but for opposite signs of the couplings g i they are {0, 1, 1}. If we assume that each end of the chain obeys one set of Wilson loop invariants we arrive to the number of end modes (1, 2, 1) in agreement with our numerical calculations. The end modes need to be symmetrically distributed in the gaps because the Hamiltonian obeys a chiral-inversion symmetry CĤ k (g, g)C −1 = −Ĥ −k (g, g) with C = σ x ⊗ σ y , which makes the spectrum particle-hole symmetric.
Global Berry phase.-The global Berry phase is defined as [62] where |ψ l,R[L] denotes the normalized l-th right [left] eigenstate of the NH bulk Hamiltonian (2). The global Berry phase is quantized to be an integer but it depends on the gauge choice, and therefore it cannot describe any observables. The discontinuities of the global Berry phase in a fixed gauge can correspond to topological phase transitions in NH systems [62], but the utility of this property is limited by the fact that the discontinuities may also appear due to the breakdown of the gauge convention. Despite these problems, the global Berry phase has been utilized in many papers [17,31,34], so we discuss it also in the context of our work. We use a fixed gauge where the first component of the wave function is real and positive. With this (optimally chosen) gauge convention, the numerical value of the global Berry phase (20) in the parameter space (g 2 , θ) with g 1 = |g 2 | leads to a good agreement with the numerically calculated number of end states (see Fig. 5 (a)). The value of global Berry phase W = 2 corresponds to two topological states at each end, and W = 1 to one pair of end states. Interestingly, the global Berry phase is not only able to identify the points in the phase diagram where the gap between different bands closes but it changes also at the topological phase boundaries corresponding to the enforced exceptional points. This happens accidentally because our gauge convention breaks down at the values of the parameters where the first component of the wave function vanishes. In the vicinity of the phase boundaries the numerically calculated W takes intermediate non-quantized values due to the limited accuracy of the numerical integration.
We further go beyond the symmetric case and study the situation when g 1 = g 2 . The results for W are shown in Fig 5 (b). In this case the results of global Berry phase calculations are not in agreement with numerical determination of number of end states. Therefore, we conclude that even with optimally chosen gauge convention the global Berry phase gives only partial description for the number of topological end states.