Interaction and motion of vortices in nonequilibrium quantum fluids

We study numerically the motion of vortices in nonequilibrium Bose–Einstein condensates, that are described by a generalized Gross–Pitaevskii equation. We analyze how the vortex properties are modified when moving away under deviation from equilibrium. We find that far from equilibrium, the radial component dominates over the azimuthal one in the distribution of vortex currents at large distances from the vortex core. The modification of the current pattern has a strong effect on the vortex–antivortex interaction energy, that can become entirely repulsive. The vortex trajectories are also strongly affected by the driving and dissipation. Self acceleration of vortices is observed in the strong nonequilibrium case.


Introduction
Topological defects play a crucial role in determining the properties of ordered phases of matter [1]. In twodimensional superfluids, the normal to superfluid phase transition is governed by the Berezinskii-Kosterlitz-Thouless (BKT) phase transition, that consists of the binding of vortex-antivortex pairs and leads to an abrupt jump in the superfluid stiffness at the critical temperature. It has been observed experimentally in experiments on Helium films [2], as well as in Bose-Einstein condensates of ultracold atomic gases [3].
Progress in the fabrication of semiconductor heterostructures has led to the appearance of another platform to study two-dimensional quantum fluids: microcavity polaritons [4]. These are hybrid light-matter quasiparticles that result from the strong coupling between a cavity photon and a quantum well exciton. Thanks to the excitonic component, the microcavity polaritons show significant interactions, where their light component gives them a light mass and allows for the straightforward optical detection through the light that leaks through the cavity mirrors. Besides the convenience for the observation, the leakiness of the cavity mirrors affects the physics because of the finite polariton life time (of the order of 100 ps at most). In order to obtain a steady state polariton fluid, some external driving-in coherent or incoherent form-is needed in order to compensate for the losses.
The first experimental detection of vortices in polariton condensates was reported by Lagoudakis et al in an experiment under incoherent excitation [5]. Contrary to their equilibrium counterparts vortices may be present in the steady state due to an interplay between disorder, driving and dissipation that lifts the constraint on time reversal invariance. Time resolved images of vortices under pulsed excitation were reported in [6], showing evidence for the annihilation of vortex-antivortex pairs. The scenario of pulsed excitation is that of the Kibble-Zurek mechanism, with a fast ramp through the phase transition. So far, no experimental evidence of Kibble-Zurek scaling was found, but it has attracted theoretical attention in several works [7,8], and stimulated the theoretical study of the subsequent phase ordering kinetics that is governed by vortex-antivortex annihilation [9]. The motion of vortex pairs was analyzed numerically in [10] in order to interpret experimental measurements on polariton condensation.
Also under resonant excitation, several works addressed polariton condensates with vorticity. Hydrodynamic nucleation of vortices was observed in [11,12]; injection of vortices in a polariton parametric oscillator is discussed in [13,14], where the dissociation of giant vortices into elementary ones was observed. Vortex-antivortex merging was reported in [15].
Under incoherent excitation, the phase of the polariton field is free and a U(1) symmetry breaking phase transition, that is the nonequilibrium analog of Bose-Einstein condensation, occurs as a function of pumping Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
power [16]. Natural questions then involve how the driven-dissipative nature of the condensate affects the superfluid [17,18] and coherence properties. For the latter, it was shown that the phase dynamics for nonequilibrium condensates is of the Kardar-Parisi-Zhang [19] type [20][21][22][23].
In a two-dimensional system, the equilibrium transition is of the BKT type and recent theoretical works have addressed the fate of this transition out of equilibrium. It was found in numerical simulations [24] that at low density, vortex-antivortex pairs proliferate, where at higher density they bind together, resembling the equilibrium case. A renormalization group study on the other hand found qualitative differences at large scales [25,26], due to a modification of the vortex-antivortex interactions out of equilibrium. Some experimental indications for a transition between exponential and algebraic decay of correlation functions have been obtained [27,28]. A transition driven by spatio-temporal vortices in a 1D system was predicted by He et al [29].

Single vortex structure
We describe the nonresonantly excited two-dimensional polariton condensate by the generalized Gross-Pitaevskii equation that is obtained after adiabatic elimination of the exciton reservoir [30,31]   y y y g y Here m is the effective mass and the contact interaction between polaritons is characterized by the strength g.
The imaginary term in the square brackets on the right hand side describes the saturable pumping (with strength P and saturation density n s ), that compensates for the losses (γ). The physical origin of the pumping term for exciton-polariton condensates is an excitonic reservoir that is excited by a nonresonant laser. For a positive polariton-polariton interaction constant g, equation (1) can be written in the dimensionless form y y y n y y . The conservative Gross-Pitaevskii equation is recovered in the limit  c 0, but also for large ν, the role of the imaginary term is expected to become negligible. We therefore have two tuning parameters that control the nonequilibrium condition of the system. We will mainly investigate the effect of c, but will also illustrate the correspondence between increasing c and decreasing ν. We assume throughout this work that the polariton condensate is in a regime that does not suffer from a modulational instability [32].
Equation (2) is solved numerically using the approach described in [33] and successfully applied to describe different kinds of vortex dynamics in superconducting condensates (see, e. g., [34,35]). We use a uniform twodimensional grid with the step = = h h 0.5 x y . The time step h t is automatically adapted in the course of calculations. This adaptation is aimed at minimizing the number of steps in t and-at the same time-at keeping the solving procedure accurate. Typically, the step h t is~-10 5 -10 −3 depending on a specific distribution of the order parameter.
The complex Ginzburg-Landau equation (cGLE) [36] is obtained from equation (2) by expanding the denominator in the last term in square brackets up to second order in ψ. The study of topological defects in the cGLE has attracted significant attention in the context of pattern formation out of equilibrium [36]. It was shown that driving and dissipation leads to a deformation of the phase profile of the vortex, that become spiral waves. The physical reason for this deformation of the phase follows from the continuity equation. In the vortex core, the density of polaritons is suppressed, leading to less radiative polariton losses, while the pumping strength remains the same. As a consequence, the vortex core becomes a source of extra polaritons, which propagate radially.
These outward flows are visible in the left hand side panels of figure 1, where we show the distributions of the current density given by the expression These distributions were obtained by evolving equation (2) in the presence of a vortex pinning potential provided by an ad hoc requirement y = 0 on a single point of the numerical grid, starting from an initial condition with unit circulation. The requirement y = 0 in the pinning point is applied after each time step integration (unless the pinning cite is removed at a certain time moment as it is considered in section 4). Since the present work is focused on the effects related to the vortex -(anti)vortex interaction and the interaction of a vortex with its own current field rather than on the influence of a specific sample shape, the boundary conditions are taken to be of the Neumann kind to minimize finite size effects. For a nonequilibrium condensate, the use of the Dirichlet boundary conditions leads-depending on their physical realization-to the appearance of additional inwards or outwards currents which would make our analysis of the vortex behavior more intricate.
From figure 1, it is clear that a larger deviation from equilibrium (c increases from top to bottom) leads to larger outward flows. At c=0.1 the radial currents are rather negligible so that a typical pattern of circulating vortex currents is seen in figure 1(a). However, for relatively large c, although the current pattern in the close vicinity of the vortex core is still determined by the azimuthal current component (see the rhs panels of figures 1(b) and (c)), the radial current component becomes dominant at larger distances from the core (see the lhs panels of figures 1(b) and (c)). In other words, for the radial current the decrease of its magnitude with increasing distance from the core is significantly slower than that for the circulating vortex current that, for all values of c, shows the typical r 1 decay at large distances [37]. This behavior correlates with the suppression of the condensate density n in a relatively wide region around the vortex core for the larger values of c (see figure 2).

Vortex interactions
A central issue that determines the coherence properties of two-dimensional quantum fluids are the interactions between vortices of the same and opposite polarities. At equilibrium, the logarithmic attraction between vortices and antivortices is responsible for their binding below the BKT temperature, resulting in a slow, algebraic, decay of the spatial correlations and a non-vanishing superfluid fraction. The logarithmic interaction potential originates from the kinetic energy due to the azimuthal flows around the vortex core. For like vortex polarities, the repulsion is responsible for the organization of multiple vortices in Abrikosov lattices [37]. As described in the previous section, out of equilibrium, the flows out of the vortex core region contain a radial component as well. It was shown analytically by Wachtel et al that the vortex-antivortex interaction acquires a repulsive component at large distances in the limit of weak nonequilibrium [25]. In this section, we will explore the vortex-vortex and the vortex-antivortex interaction energy by direct numerical simulation for arbitrary values of the nonequilibrium parameter c.
In order to compute the interaction energy, we introduce two pinning potentials, that keep the vortices (antivortices) at a fixed position. In the next section, we will describe the subsequent motion when those pinning centers are removed, which will reveal quite complicated dynamics. For a condensate with two pinned topological defects the steady-state solution of equation (2) takes the form , is the local condensate density and w + 1 2 is the system energy. The interaction energy between the topological defects is numerically determined (up to a constant) by tracking the condensate energy (frequency w + 1 2 ) for different distances between the pinning centers located at For strong nonequilibrium (large c), the radial currents dominate at large distance (see figure 1) and the current patterns for vortices and antivortices are almost the same. As a consequence, the interaction between a pair of vortices or between a vortex-antivortex pair becomes qualitatively similar. This situation is illustrated in the main panel of figure 3, where the frequency of a condensate with two interacting topological defects, w + 1 2 , is plotted as a function of the distance between them for strong nonequilibrium c=1, corresponding to the lower panels in figure 1. While the repulsion between vortices with the same sign is stronger than for opposite signs, a repulsion is observed at all distances for both polarities. We were unable to compute the interaction energies at very small distances, due to spontaneous depinning of the vortices that occurred for the considered parameters at a distance -»  figure 3) the vortex-antivortex interaction appears attractive at short distances and repulsive at large distances, in agreement with the calculations in [25]. The absence of a repulsive regime for the weaker nonequilibrium at c=0.3 could be due to the finite size of our simulation area.

Vortex pair motion
At equilibrium, the vortex interaction potential directly determines the vortex motion. Simple results are obtained for both vortex-vortex and vortex-antivortex pairs. In the former case, the pair follows a circular orbit, while in the latter case they move in parallel, analogously to an exciton in a magnetic field [38]. This dynamics can be interpreted in terms of the interaction potential and the effective magnetic field felt by a single vortex [39,40]. An equivalent explanation is obtained by considering the advection of the vortices in each other's flow field. The deviations of the interaction potential and flow fields makes us expect that these orbits will be changed out of equilibrium.
In order to investigate this issue numerically, we initialize a vortex pair at a given position with the help of the pinning potential, which is switched off after the system has reached a steady state. The subsequent motion is visualized with the help of the parameter ( ) S x y , , defined as the root mean square of the rate of changes in the local density [41] The parameter S at a given point (x, y) strongly increases when this point is passed by a moving vortex core. Figure 4 illustrates the effect of vortex-(anti)vortex interaction on vortex motion in the case of a relatively weak deviation from equilibrium ( n = = c 0.3, 1). Two vortices (vortex and antivortex), initially pinned at the positions indicated by white circles, start to move at t=0 when the pinning potentials are removed.  Let us start by discussing the motion in the vortex-antivortex case for short times, shown in figure 4(c). The initial motion qualitatively resembles that of the equilibrium case, but a deviation from parallel trajectories is visible. This deviation can be understood from the additional radial component to the current field of each vortex (see figure 1) that results in a divergence of the trajectories. For longer times, the trajectories bend due to the interaction of the vortices with the edges of the simulation area. The mechanism for this interaction is similar to that for the vortex-(anti)vortex interaction: the Neumann boundary condition at the edges impedes the radial flow, which leads to an increase of the kinetic energy. Panel (d) shows the time evolution over a longer time span, from which it is seen that the vortex and antivortex positions go toward a stable attractor, where the forces that they exert on each other are balanced by the forces due to the edges.
The motion in the case of like vortex polarities is represented in panels figures 4(a) and (b). The initial stages of the trajectories resemble the circular orbits of the equilibrium case, but again the radial flows distort the trajectories and they are pushed away from each other. The interaction with the edges of the sample stabilizes the orbit for the system; the deformation from a circular orbit is due to the square form of the boundaries. In other words, the attractor changes from two fixed points for V+AV case to limit cycle for V+V case.
For stronger deviations from equilibrium ( n = = c 1, 1), the effect of vortex-(anti)vortex repulsion due to radial currents becomes much more pronounced, as is evident from comparing figures 5(a) and (b) to figures 4(a) and (b). In case of like polarities, the circular motion of figure 4(a) is now replaced by straight trajectories up to the interaction with the boundary. For the opposite polarities in figure 5, the trajectories are no longer close to parallel: the repulsion caused by the radial flows clearly dominates over the azimuthal flow.
By comparing the time scales in figures 4 and 5, it appears that the vortex motion is faster by an order of magnitude in the case farther from equilibrium. Moreover, during the time intervals between vortex 'reflections' from the square boundaries the vortex motion may look, at first glance, 'inertial-like'. This would be rather surprising, because in typical quantum fluid situations, the inertial forces can be neglected. The parallel motion of a vortex-antivortex pair is simply obtained by equating the 'Coulomb' force (due to the interaction potential) to the 'Lorentz' force (due to the effective magnetic field that is seen by the vortex), without accounting for inertial forces. Inertial effects do exist, but require a rather subtle analysis for them to become apparent [38].
A closer analysis reveals that the 'inertial-like' ballistic trajectories are in our case not related to the vortex mass, but rather to a slow relaxation of the nonequilibrium currents. To some extent, this situation may resemble multivortex superconducting condensates, where 'inertial-like' vortex motion has been shown to arise from a finite relaxation time of the condensate rather than from an artificial mass, ascribed to fast moving vortices in earlier studies (see [34] and references therein). In order to appreciate the effect of the nonequilibrium currents on the vortex motion, we show in figure 6 the trajectory, the velocity and two snapshots of the flow field for a single moving vortex. In the panels (b) and (c), a clear anisotropy of the velocity field is visible. The orientation of the unidirectional anisotropy axis and the anisotropy strength correlate with the direction of the vortex motion and the absolute value of its velocity (see the inset to figure 6(a)). The faster is vortex motion, the stronger is the anisotropy of the velocity field around the vortex core. In other words, the vortex velocity is determined by the average supercurrent in the vortex-core region. This implies that if the average superflow through the vortex-core region is nonzero, this will affect the vortex: by 'dragging' it along.
The second ingredient that is needed to explain the apparently ballistic vortex motion is the slow relaxation of the nonequilibrium currents. In order to decouple this phenomenology from the conserved azimuthal motion, we performed a simulation without any vorticity in the polariton field. Radial currents were induced by applying an initial condition y = = < | x y t 0, 0 , which creates a core-like density suppression in the center of the square under consideration. The relaxation of currents after removing the aforementioned boundary condition at t=0 is illustrated in figure 7, showing that the response of the current to the removal of the 'core' is strongly retarded. The corresponding plateaus of j x (t) are especially pronounced for the larger distances x from the 'core'. Shorter plateaus for smaller x are partly masked at  < t 0 10 by non-monotonic perturbations of j x (t), related to the fast density recovery in the central area after removing the condition y = = | x y 0 . As further implied by figure 7, with increasing c not only the initial densities of radial currents increase but also their relative decrease with time becomes less pronounced.
With these two ingredients (drag and slow current relaxation), it becomes clear that the vortex motion out of equilibrium is complicated: when the vortex moves, it will feel its own radial flow pattern, that has not fully adapted yet to the new vortex position. We now see that this will lead to self-acceleration of the vortex core, an effect that is present even for a single vortex and is well known for spiral waves solutions of the cGLE [42]. It is illustrated by the inset of figure 6(a), which shows the vortex velocity as a function of time: the vortex is accelerated when it moves from the central region of the square towards its corner (except for the close vicinity of the boundary).
As illustrated in figure 8, the presence of radial currents significantly influences vortex motion even in the absence of vortex-(anti)vortex interactions. While at n = = c 1, 0.5 a depinned vortex moves with a gradually decreasing velocity towards the confinement-defined equilibrium position in the center of the square (see figure 8(a)), at n = = c 1, 1 no general trend to vortex deceleration appears and the vortex moves continuously with a pronounced acceleration on nearly straight pieces of the trajectory (see figure 8(b)). As seen from figure 8(c), on a long time scale the vortex motion may become quasi-stationary. The shape of the vortex trajectory for n = 1 and a rather large value c=1.5 (see figure 8(d)) shows that at very strong deviations from equilibrium the vortex motion can be affected not only by retarded radial currents but also by a delay in the relaxation of the circulating vortex currents in the vicinity of the moving core, which results in a change of direction of the vortex trajectory.

Conclusions and outlook
We have numerically studied the properties of the motion of vortices in nonequilibrium condensates. As it was known from studies on the cGLE, the nonequilibrium condition strongly affects the currents around the phase singularity. In particular at large distances, the flows becomes radial as opposed to the azimuthal flow in equilibrium condensates. The modified current pattern has a strong effect on the vortex-vortex and vortexantivortex interaction energies. For strong nonequilibrium, their interaction becomes repulsive at all distances that we were able to simulate with the help of the pinning potential. The modified interaction energy affects the motion of pairs of (anti)vortices. For example, for moderate nonequilibrium, the trajectories of a vortexantivortex pair are no longer parallel. For stronger nonequilibrium, the vortex motion becomes complicated, due to a self-acceleration mechanism: vortices are dragged by their slowly relaxing flow field. This leads to highly complicated ballistic-like trajectories, even for a single vortex. In the case of a pair of vortices with opposite vorticity, we did not see their anihilation, but further studies with a larger number of vortices will be required to understand this physics as a function of the deviation from equilibrium. Furthermore, the repulsive vortexantivortex interaction suggests a significant modification of the BKT phase transition for strong nonequilibrium. The results for multivortex nonequilibrium condensates and BKT physics, on which we are currently working, will be reported elsewhere.