Polygons of quantized vortices in nonlinear photonic waveguides

. In a nonlinear optical waveguide with defocusing Kerr-type nonlinearity, we discuss the existence of a type of stationary nonlinear waves with propagation-invariant density proﬁles, consisting of vortices located at the vertices of a regular polygon with or without an anti-vortex at its center. These polygons rotate around the center of the system and we provide approximate expressions for their angular velocity. We have computed the evolution of the vortex structures and discuss their stability and the fate of the instabilities that can unravel the regular polygon conﬁgurations. Such instabilities can be driven by the instability of the vortices themselves, by vortex-antivortex annihilation or by the eventual breaking of the symmetry due to the motion of the vortices.


Introduction
Optical vortices are regions of space in which the flow of the phase revolves around a certain point.Their study is of great importance in different types of systems and in particular for the dynamics of nonlinear fields.In the present work, we deal with vortices that are hosted in a continuous monochromatic laser beam propagating in a nonlinear step-index optical waveguide with a Kerr defocusing refractive index.Assuming linear polarization, the evolution of the amplitude of the light field ψ is given by a nonlinear Schrödinger equation of the adimensional form: where z is the propagation direction, ∇ 2 ⊥ is the transverse laplacian and Θ is a Heaviside step function describing the refractive index change of depth ∆n at the boundaries of the waveguide.The spatial dimensions are measured in units of the laser wavelength in the material and the beam amplitude is expressed in units of the inverse of the Kerr coefficient of the refractive index.
The ground state of the previous equation can be easily found by the method of propagation in imaginary time and consists of an almost constant plateau that rapidly decays to zero near the trap boundary.Vortex solutions contained in a constant background take the form ψ = e −iz e −ilθ ψ(r), being l a non-zero integer corresponding to the topological charge and ψ a real monotonic function.

Vortex polygons
We can define multi-vortex configurations as initial conditions for our simulations.Considering N V phase dislocations placed at (x j , y j ) for j = 1, . . .., N V , within the ground state background, the wave function is built as: where r j = √ (x − x j ) 2 + (y − y j ) 2 and θ j = arctan[(y − y j )/(x − x j )] are sets of polar coordinates centered at each of the vortices (or antivortices) and ψ 0 (r) is the ground state wave function of the system.The simplest shape preserving vortex structure is that of a regular polygon concentric with the symmetry axis of the waveguide.We concentrate here on the l = 1 case since it is well known that l > 1 vortices are unstable for the cubic nonlinearity [1].

Vortex polygon rotation
In the present setting, there are three kinds of effects that affect the motion of the vortices: the phase gradients generated by the rest of the vortices, effects related to the boundary of the trap and effects due to the non-point-like structure of the vortex cores, which therefore affect the condensate density in its surroundings.The velocity induced by a vortex a located at ⃗ r a = (x a , y a ) on a vortex b located at ⃗ r a = (x a , y a ) is given by [2]: where ⃗ r ab = ⃗ r b − ⃗ r a is perpendicular to ⃗ v ba .With the previous expression, we can find the velocity induced on a vortex of the polygon by the rest of them.For the a = 0 vortex, a simple computation yields R (0, 1).If we add the effect of a vortex of charge l C at the center of the polygon, we find the angular velocity induced by the phase gradients as: In figure 1, we show several cases of numerical results coming from the integration of Eq.1 with initial conditions provided by Eq 2. The l = 1 vortex polygons rotate counterclockwise, as the dominant effect on the vortex motion when the vortices are far from the boundary and from each other.It is obvious from the plot that, as expected, the angular velocity grows with N for a given R.

Static solutions
We will concentrate here on the interesting case of adding a central antivortex, since it introduces an opposite contribution to the overall rotation, allowing for transitions between clockwise and counterclockwise, leading also to the possibility of solutions with a static density profile.From our numerical simulations, we have observed that for polygon radius R close to the boundary, a counterclockwise rotation is obtained.On the other hand, for smaller polygon radius, clockwise rotation takes place.Therefore, for a particular value of R it is possible to have static solutions with zero angular rotation velocity.Moreover, numerical evolution shows that these structures are stable for long propagations.One of these solutions with R = 14 is depicted in Fig. 3, where we plot the squared modulus and the phase of ψ for a propagation of z = 80000.Without the central antivortex, the polygon would have made more than a hundred full turns during this evolution, but the plots in the figure are hardly distinguishable to the eye from the initial conditions.

Conclusions
We have provided visual insight on the dynamics of optical vortices in nonlinear waveguides through a number of numerical simulations.A more rigorous analysis of instabilities would be worthy, but it lies beyond the scope of the present work.It would be interesting to advance, both theoretically and/or experimentally, towards realistic implementations of the effects shown.

Figure 1 .
Figure 1.Rotation of a triangle (frames a-c), a square (frames df) and a pentagon (frames g-i).Dashed lines are included to guide the eye, marking the position of the central vortex.The arrows emphasize the direction of the polygon rotation.Columns from left to right correspond to z = 0, z = 100 and z = 200.

Figure 3 .
Figure 3.The stable and static configuration with a triangle of l = 1 vortices surrounding a central antivortex l C = −1, for a waveguide radius mR T = 56.The figure shows the squared modulus profile (a) and phase distribution (b) computed for z = 80000.The size of the plotted window is 140 140.The radius of the triangle R = 14 was chosen such that the positive and negative contributions to the angular velocity cancel out with each other.