Vortex creation without stirring in coupled ring resonators with gain and loss

We present study of the dynamics of two ring waveguide structure with space dependent coupling, linear gain and nonlinear absorption - the system that can be implemented in polariton condensates, optical waveguides, and nanocavities. We show that by turning on and off local coupling between rings one can selectively generate permanent vortex in one of the rings. We find that due to the modulation instability it is also possible to observe several complex nonlinear phenomena, including spontaneous symmetry breaking, stable inhomogeneous states with interesting structure of currents flowing between rings, generation of stable symmetric and asymmetric circular flows with various vorticities, etc. The latter can be created in pairs (for relatively narrow coupling length) or as single vortex in one of the channels, that is later alternating between channels.

We present study of the dynamics of two ring waveguide structure with space dependent coupling, linear gain and nonlinear absorption -the system that can be implemented in polariton condensates, optical waveguides, and nanocavities. We show that by turning on and off local coupling between rings one can selectively generate permanent vortex in one of the rings. We find that due to the modulation instability it is also possible to observe several complex nonlinear phenomena, including spontaneous symmetry breaking, stable inhomogeneous states with interesting structure of currents flowing between rings, generation of stable symmetric and asymmetric circular flows with various vorticities, etc. The latter can be created in pairs (for relatively narrow coupling length) or as single vortex in one of the channels, that is later alternating between channels.

I. INTRODUCTION
Coupled microrings (microdiscs, or more generally microcavities) are standard basic elements in diverse physical applications. In optics they are used for nonreciprocal devices [1], switches [2], loss control of lasing [3], and ring lasers [4,5], to mention a few. Recently the attention was also turned to more sophisticated devices which include several coupled microcavities and which are referred to as photonic molecules [6]. Coupled non-Hermitian microcavities are also used for study of chiral modes in exciton-polariton condensates [7], as well as for modeling coupled circular traps for Bose-Einstein condensates (BEC), where gain corresponds to adding atoms while nonlinear losses occurs due to inelastic two-body interactions. They can also be realized in nanoplasmonic systems [8] and slow-light optical microcavities [9].
Like in case of any coupled subsystems, the characteristics of coupling between microrings (or to an external element, for instance to a bus fiber) may strongly affect the stationary regimes as well as the dynamics supported by the system. The coupling can be modified in various ways. It depends on the geometry (i.e. on the mutual locations of the rings), on the wave-guiding characteristics of the rings which determine the field decay outside the cavities, on the medium between the cavities (it can be homogeneous or gradient; active, absorbing, or conservative), etc. Thus, it is of natural interest to understand how the characteristics of coupling affect the field distribution and dynamics inside the ring cavities. This is the question which is addressed in the present work. The main emphasis is made on the interplay between the size of the coupling domain and other spatial scales of the system (i.e. on the ring lengths and on the characteristic scales of the excited modes).
Before we move on to specific model (formulated in Sec. II) whose stationary solutions are investigated in Sec. III, the dynamical regimes in Sec. IV and Sec. V, and vortices in Sec. VI, we would like to advertise our effort to design up a nanostructure of coupled rings to explore experimentally this and similar settings and its dynamics. In Sec. VII we discuss specific experimental proposal based on the structures that were already manufactured.

II. THE MODEL
In the present study we consider a model described by two coupled nonlinear Schrödinger equations with gain and nonlinear loss (depending on applications they also can be termed Gross-Pitaevski or Ginzburg-Landau equations), which we write down in scaled dimensionless units Here ψ 1 and ψ 2 are the fields in the first and second waveguides, γ is the linear gain and Γ is the nonlinear loss, both considered constants along the waveguides, and J(x) is position depending coupling. Extended discussion of model (1) and of its applications can be found in a previous publication [10] where the rings were homogeneously coupled, i.e. where it was assumed that J(x) is constant. We also mentioned that (1) with constant coupling is analogous to the model introduced earlier in [11] where it was considered on the whole axis subject to the zero boundary conditions. In this paper we focus on expanding the study of the model through introducing coupling modulation J(x). We consider rings assuming, without loss of generality, that x ∈ [−π, π]. This implies periodic boundary conditions for both channels: ψ i (x, t) = ψ i (x + 2π, t), and the coupling function J(x) is extended only in the certain region of the rings. In particular, for numerical simulations we shall consider local Gaussian coupling in the following form where w is the width of the coupling, while J 0 characterizes the coupling strength. Our results are not sensitive to the particular shape of the wavefunction, as we have checked using supergaussian functions with very high power. An important remark about the used terminology is in order. For all applications mentioned in the Introduction, the meaning of the variable x is an angle defining a point on the circumference. The functions ψ 1,2 are rather envelopes of the field distributions than the total fields (see e.g. [12] for optical resonators and [13,14] for BEC applications). Thus, solutions for ψ 1,2 having nonzero topological charge (see [10]) may correspond to the total field distributions having phase singularities in the centers of the rings. In other words, such solutions describe vortices. Taking this into account the respective solutions are referred below as vortices. System (1) is simple, but possesses surprisingly rich and diverse set of stable states (some of them nonstationary). For the limit of very wide coupling (w π) we expect the same results dynamics as for the constant coupling (it is described in [10]). On the other side, very narrow coupling (w π) allows to approximate coupling with delta function. Most of the results found in the present study are numerical (using propagation techniques). In this context there is one important issue that we need to address before we present the outcome of our investigations. As discussed in [10], for uncoupled case (J 0 = 0) one can find stable background solutions in the form Once the rings become coupled, modulation instability occurs mostly due to the interplay between gain and nonlinear absorption. In the case of constant coupling in Ref. [10] two distinct classes of solutions have been found analytically: symmetric, characterized by ψ 1 = ψ 2 , and anti-symmetric with ψ 1 = −ψ 2 . The anti-symmetric solutions are always stable, and symmetric are usually unstable. Therefore we decided to perform numerical studies using symmetric state as initial condition. This led us to plethora of new states and attractors [10]. We found that even if the coupling is not uniform, dynamics starting from anti-symmetric states lead to antisymmetric stationary stable solution (see the examples of such attractors in Fig. 1). On the other hand, starting from a symmetric state, in some regions of parameters, does not necessarily lead to the antisymmetric stationary states, but can traverse to the more interesting attractors, like limit cycles or even chaotic states. Hence below we focus on the dynamics starting with the symmetric initial state with small perturbation imposed in the form of where the perturbation β was typically of order of 10 −2 . In our simulations we took the value of the gain parameter γ larger than loss Γ and we checked that all results are qualitatively the same, regardless of the particular values of this parameters. The results also do not depend on particular values of the amplitude of the perturbation β or perturbation wavenumber k. In the case of stronger loss (Γ > γ) results seem to be different and they are not included here. In this article we assume γ = 3, Γ = 1 for all later considerations and propagate from initial state defined with (4), unless stated otherwise. All simulations were performed via Split-Step Fourier method [15].

III. STATIONARY SOLUTIONS
When the coupling between the rings is weak (J 0 ≤ 1), starting from the initial state (4) we observed that the propagation leads to stable, stationary anti-symmetric solutions. Resulting wavefunction has form of constant background with the bulge in the coupling region, which depends on coupling function J 0 , as presented in Fig. 1. In this figure we plot the modulus of the wavefunction for unitary width (w = 1) and various coupling constant J 0 ( Fig. 1 left panel). In the right panel of Fig. 1 we fix J 0 = 1 and change w, going towards narrow distributions, to show what one can expect when the coupling has a form of Dirac delta function. We also show background level, plotting it as a black horizontal reference line, in left panel of Fig. 1.
In order to interpret the results of Fig. 1, we first notice that for antisymmetric solutions ψ 2 = −ψ 1 the coupling plays the role of the linear potential well, −J(x), thus leading to the equation Let us now rewrite this model in the hydrodynamic form, introducing the amplitude distribution ρ(x) as well as the For the stationary solution we obtain from (5): Due to the parity symmetry of the problem ψ 1 (x) = ψ 1 (−x) and taking into account the continuity of the solution (i.e. of the functions ρ(x) and v(x), and of their derivatives) we have the relations.
Now we observe that the maximal density ρ max = ρ(0), is achieved at x = 0, i.e. in the point of the minimum of the potential energy landscape −J(x). Thus, at x = 0 we obtain from (6) that ρ max = (γ − v x (0))/Γ, and thus v x (0) < 0 meaning that the energy flow is directed towards the minimum of the effective potential (what is expectable): Similarly one can analyze the point x = ±π, located at the opposite side of the ring diameter. Denoting . This is what we observe in both panels of Fig. 1. In particular, in the left panel of Fig. 1 we observe decrease of ρ with increase of J 0 /w. An interesting feature in the density distributions shown in both panels of Fig. 1 is the appearance of non-monotonic dependence of ρ(x) in the intervals x ∈ (−π, 0) and x ∈ (0, π). In particular, the minimal density is achieved in two symmetric points x = ±x m of the ring, different from x = ±π. As it follows from (6), at these points the absolute value of the velocity gradient v x (±x m ) is maximal. Since the solution separating monotonic and non-monotonic densities (in the intervals (−π, 0) and (0, π)) is characterized by ρ xx (±π) = 0 (at this solution the curvature of ρ(x) changes its sign), it is not difficult to make an estimate of the parameters for that solution in the case of sufficiently narrow coupling. Indeed in that case J(π) and v x (π) become negligibly small and one obtains from (6) and (7) µ ≈ ρ 2 (π) ≈ γ/Γ.
In following subsections we present results of our studies of the dynamics of our system with increasing coupling strength. We shall distinguish two separate regimes of very narrow (w 1) and extended (w 1) couplings. As we shall see in these limiting cases the dynamics can be of a quite different character.

IV. NARROW COUPLING DYNAMICS
In this paragraph we present results obtained in simulations with very narrow coupling function J(x), where we used Gaussian (2) of the width equal to w = 0.01. Our results are summarized in Figs. 2, 3 and 4. They are also described in the concise form below.
Depending on the value of J 0 we can distinguish several different patterns characteristic for the long (asymptotic) regular behavior. It can be stationary or periodic. In all simulations, the system reaches respective attractor state after a transient period, which depends on initial conditions and perturbation, and typically doesn't exceed t trans ≈ 50 in our arbitrary units.
For small coupling (J 0 1) asymptotically long time dynamics lead to stable anti-symmetric state shown with the black curve in Fig. 2. This state exhibits non-trivial phase structure, which is formed as a result of non-Hermiticity of our model.
For slightly larger coupling (J 0 1.5) when we propagate initial symmetric distribution for as long as it takes to reach steady state, we observe symmetry breaking (between channels) and our system approaches a stable asymmetric state, represented by blue (|ψ 1 |) and red (|ψ 2 |) curves in Fig. 2. We observe that local minima of the densities ρ 1,2 are achieved at x = π for the first ring and for intermediate points ±x m where |x m | < π for the second ring. To make a plot of the phase we excluded an x−independent component of the phase linearly growing in time. Note that hydrodynamic formulation allows us to discuss energy flow in the system. An interesting observation is that in the center of the coupling, i.e. in the vicinity x = 0, the energy flows in the two rings have opposite directions. Indeed these flows are determined by v j and in the vicinity x = 0 in the first ring with higher field density the current is directed outwards from the center (xv 1 > 0), while in the second ring it is directed towards the center (xv 2 < 0).
Finally, for larger coupling (J 0 2) our system tends to the limit cycle dynamics and we observe oscillations symmetrical with respect to the center of the coupling region. It is an interesting class of solutions, and analogous dynamics was found in the broad coupling regime (see the next paragraph). It seems to be universal and the general picture of the dynamics in this regime does not depend on the width of the coupling. Due to this universality we decided to discuss this class of solutions only for the broad coupling in the next paragraph.

V. BROAD COUPLING DYNAMICS
In the opposite limit, when the range of the coupling is comparable to the length of the ring (but not uniform yet) we also observe few classes of characteristic steady state dynamics and we classify them according to the (increasing) value of coupling strength J 0 . We present results for w = 1 and show asymptotic (in time) dynamics. In this case, since the most interesting steady state exhibits rather complex oscillations, corresponding to the limit cycle, we choose to present contour plots (see figure 3) and dynamical snapshots at various times in figures 4 and 6. Note that we always start from perturbed symmetric state given in Eq. (4), but our predictions are not sensitive to this particular choice.
As it is expected, increase of the coupling strength leads to the increase of the oscillation frequency. An interesting observation, however, is that the transient period, i.e. time necessary for establishing oscillations, also increases with J 0 . The oscillatory dynamics is symmetric with respect to both rings, with alternating field characteristics (amplitude and relative phase) exchanging after each half-period.
In the case of small coupling strength (J 0 3.5), the dynamics lead directly to the anti-symmetric state, analogous to the one presented in Figure 2. This part is identical in both cases of broad and narrow coupling.
We did not find any asymmetric stationary states in the case of broad coupling. This class of states was only present for very narrow coupling, as mentioned above. Instead we observe directly transition to the next phase described below, namely symmetrically oscillating states. However, in the broad coupling case this is not the last category identified in our research. Above this class, as we describe below, there is an extra region of the most interesting solutions, where vortices are produced in alternating manner, in one or the other channel.
Upon increasing the value of coupling strength, at approximately J 0 3.8 symmetric (with respect to the center of the Gaussian coupling function (2)) oscillations begin to develop, initially with very long period, which becomes shorter at higher coupling strength. Details of the oscillations can be identified in figures 3 and 4. In the first figure we now focus on the left panel (moving from the top to the bottom). First the contour plot of the time evolution of the amplitude and phase oscillations of ψ 1 is shown (wavefunction in the second channel is just shifted by half period), next we present the evolution of the norm in each channel (N i = |ψ i | 2 dx) as function of time and finally, we present more details, the evolution of the norm of each of the wavefunctions and the average within the full period of oscillations.
These results are complemented by the full view of the wavefunction during its half-period oscillations in the asymptotic regime in figure 4. The first and the third rows show the modulus of the wavefunctions (blue and red curves correspond to the first and second channels correspondingly), and second and fourth rows show the phase structure. We can trace the dynamics in which one of the channels develops two symmetric dips (see Fig. 4 (a1) and (b1)), that develop slowly, reach bottom (c1) and eventually flatten to make exchange with its partner from the second channel (see frames (e1)-(h1)). In parallel we show the phase structure (frames (a2) though (h2)), and perhaps the most prominent feature is very steep phase profile in frames (c2)-(d2). It happens exactly at the time when the two dip structures in one of the moduli (blue curve in (c1)) reach zero at the minimum. Notice that the solution stays symmetric all the time. This will no longer be true when we go to the higher coupling regime. We observe that after initial propagation (not shown) system develops into symmetrically oscillating state (shown in the time interval from t = 0 to t ≈ 10) and after several (typically 3 or 4, depending on initial perturbation) periods of oscillation system finally evolves into asymmetrical oscillations.
We investigated the frequency of oscillations of the periodic solutions. Results are presented in figure 5. It is a collective plot containing the results not only for our central example of γ = 3, but also two different values of this parameter (notice that we keep the value of Γ = 1). At small value of the coupling frequency grows rapidly, then a sudden jump occurs and further growth is linear. This jump is associated with yet another bifurcation (symmetry breaking), and as a matter of fact solutions marked with dashed lines no longer belong to the class described above, as it exhibits asymmetric oscillations. We will describe it in more details in the next paragraph. In our leading example of γ = 3 this phase transition occurs at coupling J 0 ≈ 4.5. In the transition region around this value of the coupling we observe that dynamics leads first to the symmetric oscillations and after the transient period of several symmetric oscillations it follows to the asymmetric oscillations, which are asymptotically stable. This transient behaviour can be identified in the right panel of figure 5, where we plot norm of the wavefunctions.
For values of the coupling above J 0 ≈ 4.5 (this particular value corresponds to Γ = 1 and γ = 3, see fig.5) phase structure becomes asymmetric with respect to the center of the coupling function, see Eq. (2). This leads to periodic, limited in time, appearance of vortex in one of the channels, which then (on the regular basis) reappears in the opposite channel, with inverted topological charge. The whole dynamics has again a form of regular oscillations and we show various phases of the evolution of the wavefunctions in figures 3 and 6. In Fig. 3 we show the time evolution of the modulus of the wavefunction and its phase on the contour plot. Note that in this regime there is single dip that appears once when a vortex is created (vortex in our case, as we mentioned above is equivalent to excitation) and again when vortex disappears, only to show up a bit later in the opposite ring. One can follow this process even closer looking at figure 6. In panel (a1) we see the dip which is just about to reach zero at around x = −1. The phase around this point is very steep (panel (b)), and vortex is created (there is a phase across the ring equal to 2π). Then the wavefunction flattens, until second dip starts to develop at around x = 1. Once second dip reaches zero, vortex disappears. Then the whole process repeats itself in the second ring in reversed order (panels (g)-(l)). This type of behavior was not observed in the case of narrow coupling, when we only have a symmetric structure. It seems that there is some distance between the edges of the coupling function necessary for this spacial symmetry to be broken.

VI. ON VORTEX CREATION WITHOUT STIRRING.
Our system is non-Hermitian and as such it does not conserve topological charge; vortices can be created during the dynamics, even if the system is rotationally invariant. It happens due to the modulational instability, and as we mentioned above in this 1D system vortex is equivalent to higher momentum excitations. Previously [10], for some values of the coupling constant we demonstrated that the system can, starting from perturbed symmetric state, arrive at the stationary states defined as Here κ is an integer number expressing the value of the topological charge. This antisymmetric state is stable against initial perturbations. There is also symmetric solution, where Ψ 1 (x, t) = Ψ 2 (x, t), but they are usually unstable and we never observed that in the asymptotic limit of our simulations. So far we could create vortices with equal topological charges in both channels. But here, when we define inhomogeneous local coupling, additionally varying in time (it is enough to switch it on for some time and then switch it off) there is even more exciting possibility of creating vortex only in one ring, accompanied by constant solution in the other. The idea is based on the results obtained for the broad coupling function, described in Sec. V. Imagine that we start from perturbed symmetric state (as we did in all the cases described in this manuscript) and we allow the system to evolve towards the regime of asymmetric oscillations. When the system is already in this regime, we will abruptly turn off the coupling between rings. Dynamics corresponding to this scenario is illustrated in figure 7. Here we present series of snapshots of both wavefunctions and their phases. In Fig. 6 in panels (a1) and (b1) the coupling is still on, and one of the wavefunctions is at the stage when vortex is created and dip in the amplitude is fully developed (reaches the bottom). In the series of panels (columns) in figure 7 the coupling between rings is off and we observe slow, independent relaxation in each of them separately. In this circumstances one of the wavefunctions preserves its vortex structure and tends to the solution with smooth modulus, defined in Eq. (10) with J 0 = 0, and the opposite ring develops wavefunction defined by Eq. (3). Turning off coupling in any other moment, while non-zero topological charge is present in the system, leads to even faster relaxation into described state. In this scenario we showed how to selectively create vortex in one ring. In the remaining part of the manuscript we present experimental proposal, where such dynamics may be investigated, and propose similar coupled ring systems that we intend to investigate in the near future.

VII. EXPERIMENTAL PROPOSAL: PLASMONIC AND METAMATERIAL EFFECTS IN ARRAYS OF NANORINGS
The effects discussed above can be studied via electromagnetic response of arrays of metallic nanorings. These form photonic/plasmonic crystals with enhanced electromagnetic response [8]. Such structures can be made by a variety of techniques, including the electron beam lithography (EBL), as well as the Shadow Nanosphere Lithography (SNL) [16,17]. This last technology provides an inexpensive route to complex periodic nanostructures, including nanorings. Figures 8 (a) and (b) show scanning electron microscopy (SEM) images of arrays of metallic nanorings (circular and c-shaped, respectively) deposited on a substrate using SNL. Such arrays can be used as a basis for corresponding arrays of coupled nanorings, which will have a very pronounced electromagnetic response depending on the physics of the inter-ring coupling and non-linear, intra-ring losses. To obtain the coupled nanoring arrays, we make two copies of an array (e.g. circular rings), as shown schematically in figure 8 (c), left panel. Each copy is an array, deposited on a transparent but lossy substrate (blue), and with the rings (orange) coated with a dielectric film (yellow). The two copies are sandwiched together, but depending on their relative offset, various configuration overlaps can be achieved. For example, one can perfectly align the rings (not shown) and therefore realize a uniform coupling (J(x) = constant). However, by shifting the two arrays relative to each other, one can achieve localized coupling as discussed above (configuration 1), or multiple localized coupling (configuration 2). In configuration 2 one has a more complicated situation, with some rings coupled only to single rings in the other array, and some to multiple rings (in both arrays). In addition, even in the configuration 1, one can achieve a pair of localized coupling regions by making horizontal shift adjustments. The local coupling can be also realized by tilting rings with respect to each other, since the strength of the coupling is proportional to the distance. One can also use inhomogeneous filling of the inter-ring space.
Experimentally realizable, more complex, configurations could be of interest, and we will study the corresponding models elsewhere.
The physics described in previous sections relies on sufficiently large non-linear absorption losses. These can be controlled by choosing highly non-linear materials for the substrates (or substrate coatings), such as organic semiconductors, e.g. acetoacetanilide [18]. Another way to control the intra-ring nonlinear absorption would be to employ lossy metals for the body of the ring, or use the c-shaped rings, as shown in figure 8 (b), with a nonlinear material coating in the opening. Such processing is possible [19]. The inter-ring coupling can be easily controlled via the dielectric coating (yellow colored layer in the left panel of figure 8 (c)). Also, if the substrates could be made conducting, the bias across the pair could control the inter-ring coupling as well, relying on the non-linearity of the current-voltage characteristics of these metal-insulator-metal (MIM) structures.

VIII. CONCLUSIONS
In this work we continued our investigation of simple ring-shaped nonlinear waveguides in the presence of the linear gain and nonlinear dissipation, adding local coupling between two channels. We have identified stationary and dynamic solutions both in case of narrow and broad coupling. These solutions represent different types of symmetry breaking, corresponding to bifurcations of fixed points (stationary solutions) and of limiting cycles (symmetric and asymmetric oscillatory solutions). Dynamic solutions, connected with vortex generation, allow us both to control channel populations and generate vortex states in the system. We propose experimental realization of our model in nanoplasmonic system.