Nonequilibrium polariton condensation in biannular optically induced traps

We report the mean field model of nonequilibrium polariton condensation in annular effective non-Hermitian potential traps, stemming from incoherent optically induced excitonic reservoirs of annular shape. We solve the linearized extended Gross-Pitaevskii equation in the approximation of two delta-function effective shell potentials for complex spectra of trapped polariton modes and calculate corresponding condensation threshold optical pumping powers. The exhaustive map of condensate quantum number transitions in the multi-dimensional space of trap parameters, including a cascade of topological charge increments, is drastically different from the single annular trap case in topology and the range of accessible condensate states.


Introduction
Exciton-polaritons as mixed light-matter quasi-particles, formed in optical microcavities in the regime of strong coupling of cavity photon and exciton modes, present a peculiar platform supporting interacting bosonic condensates and low-threshold bosonic lasing [1]. In particular, nonequilibrium bosonic polariton condensates in effective potential traps, generated with profiled nonresonant optical pumping, can be simultaneously confined and populated via stimulated scattering by incoherent excitonic reservoirs [2]. This system supports surprisingly rich physical phenomenology, including Berry phase stemming from exceptional [3] and diabolical points [4], formation of vortex lattices [5], spin bifurcations [6,7], spin ordering [8,9] and topological transitions [10] in trap lattices.
Unlike their equilibrium counterparts, such polariton condensates are not necessarily formed in the ground state of a potential trap if one of its excited states is strongly populated by the reservoir due to a better overlap with the latter [11,12]. In annular shaped traps this leads to formation of persistent counter-rotating polariton currents due to condensation in high angular momentum states [13,14]. In the nonlinear regime, where polariton interactions play a significant role, such condensates can evolve into space-time ordered phases [15] and persistent vortex states [16], whose topological charge can be controlled optically [17,18] or, as in the case of spinor rotating condensates [19], with magnetic field [20]. Overall, nonequilibrium optically trapped polariton condensates offer a versatile platform for generation of coherent light with optically controllable spatiotemporal parameters, including angular momentum and topology.
The model suitable for description of polariton condensation in annular traps depends on their size, shape, and particular applications. Elliptic traps of small radii, where condensation occurs at the ground state or at low excited states, are well described with the parabolic effective potential model [21]. In the case of larger traps the step-like effective potential model is reproduces the condensate angular momentum increasing with the trap size via a cascade of transitions from confined modes to continuum and to reproduce the superlinear dependence of the condensation threshold on the trap radius [15]. Delta function shell potential model also results in a ladder of angular momentum transitions in trap size and, in addition, allows assessing dissipative coupling of adjacent traps via overlapping evanescent condensate wavefunction tails [22]. In both step-like and delta function potential models the azimuthal periodicity of condensate density is explained in terms of counter-rotating vortex degeneracy lifted by imperfections of the cavity itself or the pumping profile.
In this work we address polariton condensation in biannular traps, formed by two concentric narrow shell potentials. Such optically induced traps demonstrated possibility of condensation in the first excited radial quantum number state [13]. We apply the non-Hermitian mean field model based on the extended Gross-Pitaevskii equation to compute complex energy spectra of the system. Investigating the spectrum dependence on the density of the excitonic reservoir forming the trap, we then calculate the lasing threshold and identify the condensate parameters at this threshold.
The paper is organized as follows. The theoretical model based on matching the solutions of extended Gross-Pitaevskii equation is described in Section 2. Section 3 describes numerical solutions of the model and presents phase diagrams of condensate parameters. Finally, the results and their implications are discussed in Section 4.

Theoretical model
The nonequilibrium polariton condensate is described with the two-dimensional Gross-Pitaevskii equation (see, for example, [23]) Here Ψ( , r) and (r) are the condensate wave function and the reservoir density respectively, is the effective polariton mass, and are the condensate-reservoir interaction parameters describing repulsion and stimulated scattering from the reservoir into the condensate respectively, and Γ is the polariton decay rate.
Here is the excitonic reservoir decay rate and (r) is the spatially non-uniform pumping power.
In the vicinity of the condensation threshold pumping the condensate density may be neglected in Eq.(2) on the reservoir density. In the adiabatic approximation the reservoir density relaxation time is assumed to be short compared to the characteristic timescales of the condensate dynamics. The density of the reservoir thus instantly adjusts to the slowly evolving condensate, which allows expressing the quasi-stationary reservoir density as ( , r) = (r)/( |Ψ| 2 + ). Substituting this expression, linearized in weakly populated condensate density |Ψ| 2 ≪ / , in Eq.(1) results in local anti-Hermitian dissipative nonlinear terms in the effective Hamiltonian. However, we are primarily interested in stationary states, where both condensed and reservoir parts of the polariton system are stationary. Reservoir density in this case inherits the spatial profile of the pump and may be approximated by a superposition of two concentric delta function rings: Given the circular symmetry of the system the variables are separated in Eq.(1) with substitution Ψ( , , ) = −i /ℏ+i ( ), where and are polar coordinates, is the angular momentum, is the complex energy, and the radial wavefunction obeys the stationary dimensionless equation Here = 2Γ /ℏ is the normalized radius and ( ), in the case of a biannular trap with the radial profile formed by two delta functions, given by where = / , = /(ℏΓ), and = /Γ is the normalized pumping power. The radial wavefunction converging at = 0 and → ∞ may be piecewise defined in three regions: The wavefunction continuity requires that I ( 1 ) = ( 1 ) and . The other couple of conditions on the coefficients may be obtained by integrating Eq.(4) in infinitesimal vicinities of 1 and 2 : 2 2 The above conditions form a linear system of equations on the coefficients of the wavefunction. It results in the relation between the complex energy and the pumping power : where the matrix is explicitly given in Appendix A.

Numerical approach
The resonance condition Eq. (9) was numerically solved using two different approaches. In the first approach, complex discrete spectra of energies were computed for varying values of pumping power . This allowed qualitatively illustrating the spectrum evolution with increasing pumping power and specifying the particular state that first crosses the condensation threshold Im{ } = 0. Alternatively, mode-specific critical pumping powers were computed for a wide range of quantum numbers by solving (9) for real-valued and . The condensation threshold and the condensate parameters were studied in depth in this second approach by identifying the minimal critical pumping among all modes. The results of both approaches are presented separately in the following subsections.

Complex spectrum computation
Here we illustrate the behavior of polariton complex spectra in biannular optically induced traps. Polariton complex spectrum in the presence of varying complex potential is discrete and may be labelled with integer non-negative angular number and radial number . Graphically represented evolution of the spectrum with increasing pumping power illustrates the competition of modes and the mechanism of polariton nonequilibrium condensation. Figure (1) shows the evolution of the complex spectrum for = 0, 1 = 1.7, = 2 / 1 = 1.1, and = 1 with increasing pumping power. Note that each value of corresponds to a series of spectral points on the complex plane. While the real part of the complex energy is monotonously increasing with pumping power, the imaginary part, corresponding to the growth rate (Im{ } > 0) or decay rate (Im{ } < 0), reaches its maximal value at a certain pumping power and decreases at further growing power. If this maximal value is positive, the mode is characterized with a critical pumping power, where Im{ } = 0 and may form the condensate if this critical power is minimal among all modes governed by the quantum numbers and . These critical powers were numerically computed and are discussed in the following subsection.
For pumping powers above the condensation threshold multiple modes can simultaneously have positive imaginary parts of complex energies, which results in mode competition among such modes. Although this regime is outside of the linearized model validity range, one may still expect the mode with the slowest effective decay rate (corrected for stimulated pumping by the reservoir) to form a stable condensate state if interactions weakly affect mode stability. If, however, interactions play a significant role, multiple competing modes can dynamically mix, resulting in complex nonlinear behaviour.

Critical pumping powers and quantum number maps
In this part, we present the numerical results obtained by solving the equation (9) for real-valued and using the Levenberg-Marquardt algorithm. We numerically computed the critical pumping values corresponding to various angular numbers in a range of trap radius values 1 for fixed values radii ratio and = / . Recovering the radial wavefunction ( ) with extracted coefficients of the piecewise definition (6), we then attributed these states with radial quantum numbers = 0, 1, 2, .... Finally, we computed the condensation threshold pumping power and the corresponding condensate quantum numbers and by minimizing critical pumping power values over all states in a region of the parameter space set by 1 and for fixed values of . Figure 2 shows the critical pumping dependence on the inner circle radius 1 for fixed values of = 1.1 and = 1. The red dot indicates crossing of two graphs corresponding to = 0 and = 1 ( = 0 for both) and illustrates switching between two condensate modes at 1 ≈ 1.1. For smaller traps the condensation threshold corresponds to the lowest mode ( = 0, = 0), while in larger traps the condensate forms at the mode ( = 0, = 1). The condensate density spatial distribution |Ψ| 2 of both modes is shown in the insets. Further condensate quantum number switching is visible for larger traps at 1 ≈ 1.6. The case = 1.1, where the two circle radii are close, is qualitatively similar to the single annular trap, where the condensate only forms at the ground radial state = 0 [15,22]. This is illustrated in Figure 3 showing the distribution of the condensate quantum numbers on the parameter plane ( 1 , ) for = 1.1. Similarly to the single annular trap case, switching cascade leads to the angular momentum increasing with the trap radius 1 and . Insets illustrate the spatial density distribution of the corresponding condensate wavefunctions at selected parameter plane points.
For higher values of , the condensate quantum number distribution map significantly changes. Figure 4 shows the condensate quantum number distribution in the case = 2. Most importantly, condensation at states with nonzero radial numbers ≠ 0 is possible in certain domains of the parameter plane, in agreement with the experimental findings of Ref. [13], which have not been reproduced in theory. At the same time, monotonous growth of the azimuthal quantum number with 1 and is replaced with a non-monotonous cascade of transitions, where both quantum numbers and change by varying increments. One may also note that the distance between neighbouring transitions on the parameter plane generally increases. Finally, in addition to continuous lines separating domains of condensate quantum numbers triple points lying on the edges of three domains emerge in drastic contrast to the single annular trap case. Nontrivial behaviour of mode switching in this case if underlined by the fact that eigenstate wavefunction antinodes do not necessarily overlap with the two potential peaks. This is illustrated in Fig. 5, showing the real parts of the trapped polariton mode energies and corresponding probability density radial profiles.
Finally, the average distance between neighbouring azimuthal nodes of polariton condensate density spatial distribution was numerically estimated in the same region of parameter space. Qualitatively, this parameter is approximately proportional to the ratio of the trap circumference 2 1 and the number of azimuthal nodes 2 . The 1 / ratio distribution is presented in Fig  6. Interestingly, although in each domain of quantum number distribution this value is linearly increasing with the trap radius 1 , at larger scales including many transitions the internodal distance is decreasing with the trap size in qualitative agreement with results of Refs. [12,13].

Discussion
Polariton condensation was modelled in biannular optically induced complex potential traps within the framework of mean-field generalized Gross-Pitaevskii equation. The model of biannular delta-function complex potential, emerging due to excitonic density generated by external nonresonant pumping, allowed quantitative in-depth analysis of polariton nonequilibrium condensation. A drastic qualitative modification of polariton condensation picture in comparison to the single annular trap case was discovered.
This model allowed us to reproduce for the first time and to give qualitative interpretation to formation of polariton condensates in excited radial quantum number states of the trap, experimentally observed in Ref. [13]. The developed numerical method allowed constructing phase diagrams of polariton condensate quantum numbers in arbitrary regions of the multidimensional space of optical trap parameters. This paves the pave towards harnessing the radial degree of freedom of polariton laser emission via all-optical control of polariton condensate angular and radial quantum numbers.
One of the most intriguing predictions of the model is the modification of the condensate quantum number distribution map topology with increasing ratio of the biannular trap radii. It supplements the layered cascade of incremental quantum number transitions with triple points in angular quantum number and quadrupule points if the radial quantum number is taken into account. Topological characteristics of such points and the possibility of exceptional point emergence in this system remain unrevealed and may be studied within the developed model framework.
It should be noted that the discussed phenomena rely on both Hermitian and anti-Hermitian terms of the effective Hamiltonian, stemming from exciton-exciton repulsive interactions and bosonic stimulated particle exchange between the reservoir and the condensate respectively. It is illustrated by qualitative difference of condensate mode switching patterns on the parameter , Fig. 4. Distribution of quantum numbers (shown with colour) and in the parameter space, formed by the trap inner radius 1 and the interaction parameter ratio = / in the case of a higher ratio of radii = 2 / 1 = 2. Spatial condensate density distribution is shown with round insets for selected points on the plane. quantifying the ratio of the two types of condensate-reservoir interaction. Being an interacting non-Hermitian system, nonresonantly pumped nonequilibrium exciton-polariton condensates present an opportunity to vary this parameter by independently controlling exciton-exciton repulsion and stimulated scattering rates with the Hopfield coefficient on the one hand, and the energy dispersion anti-crossing shape on the other. Fig. 6. Numerically computed approximation for the distance between neighbouring azimuthal nodes of condensate density = 1 / in the same region of parameter space as in Fig. 4 ( = 2).

Appendix A
The general view of the resulting matrix on the spectrum can be given as