Optical vortices in waveguides with discrete and continuous rotational symmetry

Coherent vortex structures are fascinating physical objects that are widespread in nature: from large scale atmospheric phenomena, such as tornadoes and the Great Red Spot of Jupiter to microscopic size topological defects in quantum physics and optics. Unlike classical vortex dynamics in fluids, optical vortices feature new interesting properties. For instance, novel discrete optical vortices can be generated in photonic lattices, leading to new physics. In nonlinear optical media, vortices can be treated as solitons with nontrivial characteristics currently studied under the emerging field of topological photonics. Parallel to theoretical advances, new areas of the engineering applications based on light vortices have emerged. Examples include the possibility of carrying information coded in the vortex orbital angular momentum, understood as a spatial-division-multiplexing scheme, to the creation of optical tweezers for efficient manipulation of small objects. This report presents an overview highlighting some of the recent advances in the field of optical vortices with special attention on discrete vortex systems and related numerical methods for modeling propagation in multi-core fibers.


Introduction
Coherent vortex structures are often associated with hydrodynamics and fluid motion. In a fluid vortex, the intensity distribution forms a hole corresponding to the singularity of the velocity field. In a similar manner, electromagnetic waves can create phase singularities embedded in the continuous field with zeroes of intensity distribution at such special points (pivot points). Light intensity plays a role similar to density in hydrodynamics, and the gradient of the optical phase is akin the fluid velocity. Similarity of the mathematical models makes possible the transfer of knowledge and ideas between the two fields.
Optical vortices present an interesting and growing area of research that combines fundamental theoretical and mathematical science and emerging applications (see, e.g., recent reviews ( [1][2][3][4][5][6][7][8][9][10]). Vorticity is characterized by an integer winding number (topological charge), which is *Correspondence: pryamikov@mail.ru 1 Laboratory of Spectroscopy, Prokhorov General Physics Institute of RAS, Moscow, Russia Full list of author information is available at the end of the article related to a total change of the phase during the round trip(s) over the vortex pivot point. Indeed, optical vortices can be seen as a physical realization of mathematical concepts such as topological defects, the accumulation of geometrical phases, and the presence of a phase singularity in a complex field (see, e.g., [1][2][3][4][5][6]11] and references therein). Vortex structures are realizations in the fast growing field of the topological photonics. All forms of vortices in linear and nonlinear fields obey the conservation of the topological charge S, which is defined as the total change of the phase along a closed curve surrounding the pivotal point of the vortex divided by 2π. The special behavior of the phase and amplitude near the vortex pivot point results in the circular flow of energy in optical vortices ( [12][13][14][15][16]) and this property allows optical vortices to carry orbital angular momentum. On the applications side, light vortices are of interest for various applications, such as in creating optical traps [17][18][19], modes of information transmission [7,9,20], astrophysics [21,22], microscopy [10,23] and laser micromachining [24,25]. The inherent robustness of topological properties provide for extra stability of vortex structures.
In nonlinear media, optical vortices are treated as (topological) vortex solitary waves [2,11,26]. Such topologically stable pulses also can be used as information carriers [9,27,28]. We should note that while continuum vortex solitons are highly sensitive to azimuthal instability, stabilization and arresting this instability can be achieved by different means [29][30][31]. Alternatively, an interesting possibility to create stable structures, is to generate them in optically induced photonic lattices. This unique scenario is clearly absent in a hydrodynamics setting; one that enriches the field of optical vortex dynamics. Examples of optical lattices include photonic crystals, arrays of waveguides, or those created by interference of optical beams One particularly special applications of optical vortices is its use as an information carrier. Indeed, light vortices seen as carriers of orbital angular momentum can support spatial-division multiplexing, which is an actively studied topic in high-speed optical communications addressing the continuously increased demand for capacity in global communications systems. Using spatial degrees of freedom in fiber-optical communications systems is a topic of active research and new concepts and ideas are in high demand in this field. With this motivation, multimode [32,33] or multicore fibers (MCFs) [34,35] may be used for such spatial-division multiplexing in fiber-optic transmission links. Apart from the telecommunication applications, MCF is an interesting physical system that features spatial discreteness, or a combination of discrete (space) and continuous (time) properties, as well as nonlinear effects. This photonic structure facilitates the propagation of spatiotemporal optical solitons [36], and the corresponding spatiotemporal optical vortices which were studied theoretically in [37] and experimentally in [38,39]. Beyond these early studies, distinct nonlinear effects present in MCF have been studied in [40][41][42][43][44][45] with results indicating fascinating new features when compared with both infinite photonic lattices or single core fibers. Mathematical models describing light propagation in MCFs are centered on the nonlinear Schrödinger equation [46][47][48] or the nonlinear complex Ginzburg-Landau equation [49][50][51][52]. The latter non Hamiltonian model describes an active MCF media with loss and gain effects affecting light propagation.
This review will be built around the following basic mathematical model, to be considered in its most generic form as: How this model arises is shown in some detail in "Optical light bullets in discrete systems" section. It includes 2d (n,m) or 1d (n=m) discrete configurations in space. We will examine it both for infinite lattices and in particular configurations with a finite number of sites more realistic for practical applications. Observe that depending on the importance of the temporal dynamics, governed by the third term (group velocity dispersion) in Eq. 1, this system describes 1d (not considered here), 2d and 2+1 dimensional nonlinear systems. We assume a Kerr-type nonlinearity as the last term in the l.h.s. of Eq. 1 indicates. The third term in Eq. 1 is due to the dispersive properties of the medium. For instance, in multi-core fibers, this term is responsible for the temporal evolution of the field and its corresponding pulse formation. Finally, the second term in Eq. 1 -(C U) nm represents the linear coupling features at site (n, m). It is clear that there are many coupling configurations; here we limit the analysis on two cases of interest: the uniform square and the hexagonal geometries, for which the coupling operators are respectively: hexagon nm = c(U n−1,m +U n+1,m +U n,m−1 +U n,m+1 +(2) U n,m+1 + U n+1,m+1 + U n+1,m−1 ) As stated above, this general model can be applied to both purely spatial discrete systems (no time dependence) and discrete-continuous scenarios. Practical waveguide arrays with a relatively small number of cores are evidently sensitive to boundary conditions, thus we will discuss differences between infinite and finite lattice dynamics. A principal objective we target relates to recent advances in both discrete and continuous optical vortices. Published research on optical vortices is already immense and the area is expanding fast. With this in mind, this paper is not meant to be a comprehensive overview of the field; instead, we rather focus on several topics which are close to our research interests with a special emphasis on discrete vortex structures and the underlying numerical modeling of complex light fields in MCFs. We begin with discussion on optical vortices for infinite 2D lattices in "Basics of discrete vortices in nonlinear optical lattices" section. "Discrete vortices in cylindrically symmetric nonlinear systems" section examines finite-element discrete systems with a special focus on MCFs. "Optical light bullets in discrete systems" section deals with optical light bullets in continuous-discrete systems. In "OAM modes and the poynting vector singularities in optical fibers" section, we discuss optical vortices in linear and nonlinear

Basics of discrete vortices in nonlinear optical lattices
We now introduce basic facts concerning rotating discrete localised structures in infinite lattices. Different types of discrete lattice vortex solitons have been theoretically predicted in [53][54][55][56] and experimentally demonstrated in [57,58]. Further advances in photonic technologies open the possibility of discrete vortex generation in photonic crystal fibers and other novel photonic lattices [59,60]. Some of the unique properties of discrete vortex solitons can be associated with soliton clusters, which feature interesting mobility properties and rotating propagation [41,42,59,[61][62][63]. Particular consideration is now given to discrete vortices in square and non-square lattices and in twisted multi-core arrays.

Vortices in 2d square lattices
The master equation describing optical field evolution in a 2d square lattice with Kerr nonlinearity reads: +U n,m+1 − 4U n,m ) + |U n,m | 2 U n,m = 0 derived from the corresponding Hamiltonian The system also conserves total power (also called in the literature the particle number) N = n,m |U n,m | 2 . Note that Eq. 3 has a family of stable two-dimensional discrete solitons (see, e.g. [64][65][66]). These classical discrete solitons do not have angular momentum.
While exact vortex solutions are not known, for the case of a weak coupling (c small) one can demonstrate the existence and stability of solutions to (1) using a mathematical technique known as the anticontinuum limit [61]. It can be easily seen that for the case c = 0, each lattice site is an independent entity with an exact solution U The building blocks creating vortices relies in defining a closed contour S of sites for which A = 1, while for the rest of the sites A = 0. From this, extending to small nonzero values of the coupling constant c, the existence of vortex solutions comes from the Implicit Function Theorem (IFT), which is the central piece the bifurcation theory. More precisely, for a simple closed discrete contour S on the plane (n, m) a discrete vortex is a localized stationary solution of (3) which in the limit c = 0 approaches (5). If we label the phases in counterclockwise order of a c = 0 vortex with N latice sites in S as θ n,m , (n, m) ∈ S, θ 1 , θ 2 , ...., θ N+1 = θ 1 , the condition for persistence for c > 0 small, based on the IFT is the "zero power flow" at every site. To illustrate by way of examples, for small c, vortices of charge 1 with energy localized in 4 sites (think of a square), which in the limit c → 0 have phase values θ 1 = 0, θ 2 = π 2 , θ 3 = π, θ 4 = 3π 2 . A second example for 8 sites is a one parameter family: with charge = 1 for θ = π 4 , charge = 2 for θ = π 2 and charge = 3 for θ = 3π 4 . Beyond the existence proof, the results summarized in [61] for vortex solutions with weak coupling show that those with charge one are stable while vortices with charge two and three are unstable.

Vortices in non-square lattices
Given that technology can build optical lattices with different geometries, it is not surprising then to see efforts to study vortices beyond square arrays. One realization considered in several studies is of a hexagonal array. As it relates to vortex dynamics, Terhalle and collaborators [67] experimentally showed stable vortices of charge 2 in hexagonal lattices in a photorefractive crystal. These lattices were formed by using a spatial light modulator which converts a beam into three interfering planes. The resulting interference pattern is that of a two-dimensional hexagonal photonic lattice "imprinted" in the photorefractive crystal. The experimental result, a rather surprising one, found in the case of self-focusing nonlienarity the single-charge vortex solitons are unstable. Instead, the double-charge vortices are stable and this is a consequence of the intersite power exchange in the vortex soliton. These features are reversed in the case of de-focusing nonlinearity -single-charge vortices become stable while double-charge vortex solitons are unstable.
It is important to point out that in some experimental settings the lattice is "drawn" in, for example, a photorefractive material. In such cases optical vortices can be described by continuum models, like the nonlinear Schrödinger equation in the presence of a potential representing the lattice. With this in mind, an important mathematical question is to assess the validity of discrete models that emerge from a continuous one. The master model for non-square spatial-only lattices is given by Eqs.
(1) and (2). these models and their solutions are discussed in more detail in a more complex spatio-temporal context in the "Optical light bullets in discrete systems" section.

Twisted multi-core arrays
The topological nature of discrete vortices is defined by having a site where the field is zero ("dark" sites"). This can be linked to the tunneling phenomena. Recently there has been much interest in analysis of a imposing a twist to a multicore array along the direction of propagation. In terms of the mathematical model, the twist adds a Pierels phase to the coupling constant amongst nearest neighbors [68]. For twisted arrays, recent works have shown what is interpreted as an optical Aharonov-Bohm suppression of tunneling, both in the linear [69] and the nonlinear [70,71] regime. Since the twist creates an orbital angular momentum (OAM), twisted multi-core arrays presents opportunities to explore novel dynamics [72] and find nonlinear modes including stable discrete light bullets [73] and vortices carrying an additional degree of freedom.
In the simplest case of a circular array that is N twisted in the direction of propagation the master model reads: Here the phase φ is proportional to the twist rate. To illustrate the topological nature of the twist, consider the N = 4 array with φ = π 4 in the linear case γ = 0. An exact solution reads indicating the third core remains dark throughout the propagation. The persistence of this feature in the nonlinear regime was shown in [70] for 4 cores and in [71] for 6 cores for the case φ = π 6 . Twisted multi-core arrays can be created in multi-core fibers (MCFs) and will be considered in more details below. The axial twist is created during the MCF draw process. Twisted MCFs are promising in optical fiber sensing applications. For instance, when the outer cores are twisted forming a helix around the central core, the core in the center can act as a reference along the length of the fiber. The outer cores are sensitive to bend and twist. Namely, a fiber twist strains all outer cores in the same way, whereas bending leads to different strains in each outer core allowing for three-dimensional measurements.

Discrete vortices in cylindrically symmetric nonlinear systems
Now we discuss in more detail discrete vortices assuming cylindrical symmetry of the undelrying mathematical models. The properties of vortices in the axial symmetric discrete optical systems differ dramatically from their counterparts in the uniform axial symmetric optical systems. With respect to the direction of the discretization, there are two main configurations of discrete axial symmetric optical systems for which the discretization can break the translation symmetry in a radial or azimuthal directions. The third possibility is with simultaneous discretization in both directions. These axial symmetric optical systems with different refractive index profiles in a radial and azimuth directions offer new opportunities for the creation, existence, stability properties and management of optical vortices.
The most intensively studied axially symmetric optical systems discretized along the radial directions are Bessel optical lattices [74][75][76][77][78][79][80][81][82][83] with a refractive index profile along the radial direction in a form of zero order or higher order Bessel functions R(η, ξ , ) = J 2 n (2b lin ) 1/2 (η 2 + ξ 2 ) 1/2 , where b lin is the radii of the lattice rings and n is the order of Bessel function. These profiles can be imprinted optically with non diffracting zero or higher order Bessel beams which preserve their transverse energy distribution during the propagation through the linear optical material. The simplest Bessel lattice imprinted with the zero order Bessel beam has central guiding core surrounded by multiple concentric guiding rings while the lattices imprinted with higher order Bessel beams have no central core but only multiple concentric guiding rings. Similarly, the different refractive index profiles along the radial direction can be imprinted with other non diffracting beams like parabolic beams and Mathieu beams. Periodic modulation of the non diffracting lattices breaks the translation symmetry in aazimuth directions and each of the multiple concentric guiding rings brakes into a finite number of guiding cores. Optical in azimuthally modulated Bessel photonic lattices, have received much attention [81][82][83].
An example of axial symmetric discrete optical system where the translation symmetry is broken in the azimuthal direction, is a multicore fiber with cores concentrically arranged in the form of a ring [40,[84][85][86][87][88]. Inclusion of the central core breaks also the translation symmetry in the radial direction. The inter-core coupling in multicore system is a consequence of the energy exchange between the cores via the evanescent-fields of guided light modes, whose extent could be set by the separation between cores, number of cores and index profiles of the cladding and cores [87].
Multicore fibers recently attracted a great deal of interest in the field of high-capacity communication systems Page 5 of 28 [32,34,89,90]. However, their applications can go well beyond optical communications due their ability to provide coherent transport of light that might be important for nonlinear devices [40,[84][85][86][87][88], coherent beam combining and new concepts and techniques in the field of high-power fibre lasers [91][92][93][94]. The light vortices -ring shaped modes guided by the MCF with linearly coupled circular array of cores -are shown to be promising candidates for these purposes. Light vortices in such system are dynamically stable structures characterized by invariant orbital angular momentum and energy [5,7,[95][96][97].
They can carry huge power over long distances without significant losses by properly arranging the system geometry and structural parameters of the cores and cladding. Namely, the inter-core coupling in a multicore system is a consequence of the energy exchange between the cores via the evanescent-fields of guided light modes, whose extent could be set by the separation between cores, number and index profile of the cores and cladding [87]. The MCF system studied in [86,87] is with the index profiles of the periphery cores fixed and identical, while the number of cores and their separations are variable. An alternative approach is to engineer the refractive index profile and geometry of the cores. The latter offers possibilities for breaking the axial symmetry and favoring propagation through a chosen core or cores. These can be relevant in designing long fibre lasers and other applications that crucially depend on the ability to provide coherent propagation of light over long distances. However, the material properties are becoming notably affected by increasing Fig. 1 Schematic of an MCF with M = 5 periphery cores and the central core. Geometric parameters that dominantly determine coupling coefficients C 0 and C 1 , are the distance between the central and each of periphery cores R and the distance between two neighboring periphery cores d = 2R sin(π/M). Coupling coefficients depend on the geometry and refractive index of cores. Here, all periphery cores are assumed to be identical, while the central core can be different the light intensity imposing different conditions for light propagation, which then can modify guided mode characteristics and destroy the propagation coherence. However, theory predicts the possibility that the high power light vortices can be tied up in circular MCFs and preserve coherent propagation over long distances in the presence of inevitable material losses and nonlinearity [86,87] . Furthermore, the presence of a central core in the MCF can offer the possibility to manipulate vortices traveling across the circular periphery core array.

Mathematical models
The model equations of MCFs ( Fig. 1) derived from the Helmholtz equation [98] are based on a lowdimensional version of the discrete nonlinear complex Ginzburg-Landau equation in a weak coupling approximation, assuming significant interactions between nearest neighboring cores [40,84]: Here, U 0 and U m , (m = 1, 2, ..M) stand for the complex amplitude of the light field in the central and periphery cores, respectively, with corresponding linear propagation constants β 0 and β m . The coupling between the neighboring periphery cores is denoted by the coupling parameters C m,m−1 , C m,m+1 , between the central and periphery cores by C 0,m , while the nonlinear parameters of the central and periphery cores are γ 0 and γ m , respectively. In general all coefficients are complex with the imaginary parts accounting for loss and gain processes in the MCF system [98].
In the case of passive and lossless MCF with identical periphery cores which ensures equalities β m = β 1 , C m,m−1 = C m,m+1 = C 1 , and γ m = γ 1 and identical coupling between central core and periphery cores C 0,m = C 0 the Eq. 9 can be rewritten as: Two conserved quantities characterize light propagation through the MCF system (10): the norm (total power) Conservation of P and H can be used for control of the accuracy of numerical modeling and for analytical search of solutions in the case of extra symmetries; e.g. assuming that the field dynamics in all periphery cores is similar.

Eigenspectrum of the MCF
The eigenspectrum of the MCF can be obtained by a linear stability analysis [40,84,86,87]. After an introduction of small perturbations to the stationary solutions, the following eigenvalue matrix equation was derived: where N is the total number of cores: N = M or N = M + 1 in the case without or with a central core. In general, the matrixM is non-Hermitian, and the type of eigenvalues depends on the entries of the matrix elements. The eigenvalue problem was solved analytically in two cases, for linear MCF without and with central core and nonlinear MCF without a central core.
The eigenvalue problem for a linear MCF without central core, Eq. 13 is a Toeplitz type of problem with known eigenmode structure [40,84,86,87]. By straightforward algebraic manipulations, it is possible to obtain expressions of eigenvalues (EVs): In the presence of a central core, the first M − 1 EVs coincide with those given by Eq. 14 (k=1,...,M-1), while two remaining additional EVs are given by the following expression where indices R, I denotes real and imaginary parts. A general linear eigensolution for a propagating mode can be expressed as a superposition of eigenvectors that form an orthonormal basis [87,95]. Any initial excitation of the array propagates as a linear combination of waves with phase velocities determined by Eqs. (14) and (15) where (17) are the pairs of eigenfunctions arising from the presence of the central core and a k , b k , k = 1, .., M+1 are free complex coefficients (specified by the launching conditions).
A linearly coupled MCF without or with a central core can support an infinite number of solutions, since every superposition of its eigenmodes is a solution. In the absence of losses all EVs are real and consequently solutions are neutrally stable. In general, solutions can be classified with respect to their angular momentum: zero and nonzero momentum solutions (vortices). In terms of propagation, solutions can be periodic or quasi-periodic depending on the system parameters. If the eigenfrequencies (eigenvalues) of a solution eigenspectrum are commensurate, the solution is periodic with regular revivals of the input state along the fibre. Hence, expressions (14) and (15) offer the possibility to engineer the fibre coupling coefficients by analytically solving the inverse eigenvalue problem. Analytical solutions of inverse problems are rare and extremely efficient in general, underscoring the complexity of these problems [87].

Vortex modes
Two intrinsic features of vortex solutions, which are one type of eigen-modes in the MCF considered here, are first the topological charge S, defined as the constant phase difference between all the neighboring cores in the periphery, and second, is the conservation of angular momentum. If a MCF contains a central core and supports vortices, the central core remains dark and does not take an active role in the vortex propagation.
In [87], it is shown that a lossless linear MCF with M = 3, 4, ...16 periphery cores can support dynamically stable vortices with topological charge S = 1, 2, 3, 4. The form of discrete vortices characterized by topological charge S is All linear vortices are stable no matter if either the central core is present or not.
Vortex properties in a more general system with nonlinearity included were considered in [86]. In contrast to linear vortices, the properties of nonlinear ones are power dependent. This accounts into the narrowing of the stability region of vortices. In the case of a nonlinear vortex in MCF without a central core, the cumbersome but straightforward algebraic calculations lead to the following sufficient stability conditions [86]: I) If cos(α) ≤ 0 then stable vortices are characterized by This inequality gives the following sets of stable vortices in the low power region: Numerical modeling demonstrated excellent agreement with the analytical results. Figure 2 summarizes the results of the linear stability analysis for the MCF with and without the central core.
The general conclusion is that the nonlinear MCF without a central core can support stable propagation of discrete vortices with unlimited power for particular values of parameters S and M. These results are consistent with the findings in the Refs. [53,54,61,99,100], which predicted stable propagation of vortex structures with S = 1, 2 in different two-dimensional lattice configurations. The presence of the central core enables the coupling of the perturbations in the periphery cores giving rise to instability for some intermediate vortex powers. However, a remarkable feature is that stability of high power vortices is not affected by the nonlinearity and the presence of the central core.
The developing of vortex instability is associated with the active inclusion of the central core in the energy redistribution of the system and the destruction of the vortex phase pattern -decoherence; see Fig. 3. As a consequence of the instability, the vortex structure evolves into a breathing mode with a more complex internal structure. The most significant finding in [86] is the possibility of a stable and coherent transfer of high energy through the nonlinear MCF via vortices.
The demonstrated remarkable feature of stable vortex propagation carrying high optical power and preserving the orbital angular momentum has significant potential for many applications. These results can be used to develop a new model of orbital angular momentum multiplexing and high capacity transmission lines using MCF. The considered discrete optical vortices also have the potential for high power field trapping and transfer in MCF, which provides a method to manage and control the nonlinearity and to design new type of all optical switches [55][56][57], sources of high brightness coherent radiation or uniform light sources with tailored phase profiles, generation and management spatiotemporal vortex light bullets in MCF settings [58][59][60]97], and supercontinuum generation [42].

Optical light bullets in discrete systems
Structuring of light in space and time is a fascinating area of research and technology. In space, light can be localized and configured by using waveguides that are formed by appropriate variations of the refractive index. Nonlinear optical properties give another practical possibility to localize and control light both in space and time. The combination of these two features lead to a rich variety of interconnected methods to structure and manipulate spatial and temporal properties of light. In particular, advances in fiber optics technology over the past 30 years, has had an impact not only in numerous highly important applications, but it has also provided a laboratory to display nonlinear phenomena such as modulation instability, soliton formation and interactions, supercontinuum generation, parametric amplification, optical wave turbulence and many others. We should be remiss not to refer to the early efforts creating localized wave-packets by balancing diffraction with nonlinear Kerr self-focusing. The formation of a spatial soliton opened a broad range of applications. Advancing this soliton theory and its instability in the 1+2 dimensional nonlinear Schrodinger equation, i ∂ ∂z U(x, y, z) + T U + |U| 2σ U = 0, proved that for the Kerr nonlinearity (σ = 1), a 2-d spatial soliton is unstable in that it either collapses (for powers above critical) or diffracts. In this equation, U represents the envelope of the electric field; T stands for the spatial variables (x, y) An alternative and practical mechanism of light wave collapse stabilization and the formation of localized solutions (light bullets) is to consider a discrete spatial photonic structure. In particular, if instead of a continuum field we have a discrete model with the corresponding discrete Laplacian. Our earlier pioneering work [36,101] on such a model, demonstrated in a specific rectangular geometry of waveguide arrays stable bullets propagation. As it has been the case in other instances, it was years later that this theoretical discovery was demonstrated experimentally. At the time, we witnessed technological advances in photonic crystals and multi-core fibers. While these were mostly driven by major challenges in optical communication. in finding methods and techniques capable of offering transmission capacity above the limitations of the single mode-fiber communication channels [102,103]. A multi-core fiber allows one to implement spatial division multiplexing, enabling a scale-up in transmission capacity per-fiber. Recent experimental demonstration of light bullets (LB) (predicted in [104]) in an array of optical waveguides [105,106] paves the way for broad applications of light bullets in relatively low cost MCF. As the technology of multicore fibers continues to advance, so are the possible applications of such arrays. Note, that spatial de-multiplexing is also important in emerging applications of multi-core fiber is in the field of high-power fiber lasers [107]. Here, the use of multi-core fiber allows one to split the total high power into channels with power below any undesirable nonlinear effects. In other words, laser beams in each core may be transported safely being below the threshold of the detrimental nonlinear effects while the total coherently combined power can be high. In this respect, our work represents an important contribution to this application in that it studies for the first time the stability properties of LB under more general coupling schemes. Finally, this new multi-core fiber technology opens up new perspectives for the fascinating research on light bullets (see e.g. [108][109][110][111][112][113][114][115] and references therein) and can be a natural laboratory to study fundamental phenomena such as nonlinear Anderson localization, light filaments and vortices and optical rogue wave formation to name some. The MCF is a specific realization of fiber arrays with flexible mutual arrangement of cores. It is important to understand how the mutual arrangement of fibers will affect the existence the LB and their stability, which is the subject of the present section. We first present a generic analytical description of localized bullets and complete analysis of it stability. The numerical calculations support the analytical results and demonstrates the range of applicability of the analytical approach. Finally, the formation of the stable light bullets from the pulse lunched in one fiber core has been demonstrated. Given the importance of LB in nonlinear science and applications mentioned, further explorations should be expected. With this in mind, the main objective of this section, is to analyze the optical bullet features in continuous-discrete optical media including their stability under the most general coupling schemes. Our focus is on the generic models that may be relevant to various applications. This theoretical work we believe, will provide a framework in the design of multicore elements aiming at optimizing desirable and specific applications such as routing, switching and coherent beam combining.
For arrays of waveguides where light propagates mainly in the central (core) region of the individual elements and for which we recall, transverse exchange of energy is due to tails of the field overlapping neighbor waveguides, the field is well approximated by a superposition E(x, y, z, t) = nm U nm (z, t)F(x − x m , y − y m )e i(kz−ωt) + c.c where we assume each waveguide is identical and supports a single mode F. The center of each waveguide is at location (x m , y m ). In this case the equations shown in the Introduction (eq. 1), in dimensionless form, describing the propagation of the envelopes propagating in the fiber at site (n,m) read: where (C U) nm represents the linear coupling functional form at site (n, m) obtained by integrating in (x, y) the resulting overlap of the mode centered in the (n, m) core with the tails of neighbor modes. As in previous sections, two cases of interest are (see Fig. 4): the uniform square and the hexagonal geometries, for which the respective This system has two known conserved quantities; the Hamiltonian H = nm (N(U; U) nm + | ∂U nm ∂t | 2 − |U nm | 4 )dt and the total power P = nm |U nm | 2 dt. The waveform of the discrete-continuous light bullets (that are extrema of the Hamiltonian under fixed total power P) U nm (z, t) = A nm (t)e iλz is given by the equation: Highly localized bullet solutions in some limits can be derived using an asymptotic approach. The consideration is that spatially, most of the energy is concentrated in one site, (0, 0) and small satellite pulses of decreasing amplitude propagate in subsequent layers. Mathematically, this means we seek solutions of the form U nm (z, t) = A nm (t)e iλz + cc with λ >> 1, where each envelope has an expansion of the form A nm (t) = A   functional form is the same for the central core for both square and hexagon lattices, the second line gives the functional form for the first layer in the square array, while the third gives the form for the first layer of hexagonal array It should be stated that these solutions fail to be uniformly valid beyond |t| ≥ O( √ λ). In fact at the pulse tails, all waveguides have solutions of the same order It is then natural to follow the study by comparing these analytical results with numerically computed general localized LB solutions.
Our numerical studies present solutions describing continuous discrete light bullets of Eq. (1) where we varied the number of elements in the array, which is an important consideration for practical systems for which the number of cores is finite. We find that the asymptotic functions a 0 (t), a 1 (t) coincide with numerical solutions A 00 (t) and A ±1,0 (t) = A 0,±1 (t) of Eq. (23) up to the order O( 1 λ 3/2 ). To fully assess the validity of the asymptotic analysis, we compared them with full numerical solutions for the case of a coupling constant c = 1 and for the NxN rectangular structure. As Figs. 5 and 6 show, there is good correspondence between the theory and numerical modeling for large values of λ. Figure 5 depicts the dependence of the amplitude of the solution in the central and neighboring cores on the parameter λ. Figure 6 shows comparison of the time domain structure of the light bullets for three different values of λ. As expected, the correspondence is good for a value of λ equal to 14 and fails at λ = 7, thus the limiting value for which the asymptotic state is valid is around λ = 8. Figure 7 summarizes the spatio-temporal amplitude and phase features of the light bullet.
Finally, global quantities of the LB such as power and the Hamiltonian are shown in Fig. 8. Here we observe point to the universal behavior for different number of fiber elements. Specific to the power dependence on λ, below we discuss how it relates to the stability properties of the LB. This behavior coincides with that of figure 2a in [105].
We now investigate the stability of solutions of the form U nm (z, t) = A nm (t)e iλz by linerarization U nm (z, t) = (A nm (t) + f nm + ig nm )e iλz , leading to the system of linear equations Defining inner products for square arrays: and for hexagonal arrays: In both instances and similar to the 1 + 1 + 1 case [64,116,117], the following properties of the linear operators hold: (i) < f , H + f > ≥ 0 and it is equal to 0 if f nm = 0 or f nm = A nm .
(ii) There exist some F for which < −F, H + F > < 0, < F, A > = 0. We should point out that while we only discussed two specific geometries, (i) and (ii) will be generally true for a large class of coupling schemes. Properties (i) and (ii) allow us to conclude that the existence of negative eigenvalues of the operator H − is a sufficient condition for instability of the nonlinear state. Furthermore, one can show this condition is equivalent to the Vakhitov-Kolokolov criterion on the sign of d dλ P = d dλ < A, A >, which proves instability of the left branch of solutions in Fig. 8 (left). For the highly localized solutions presented in the previous section, one finds that the power P(λ) = 2λ 1/2 + Kλ −3/2 , where the constant K depends on the coupling coefficient and the geometry. Observe that a minimum is achieved at λ c = ( 3K 2 ) 1/2 , so that stability of the localized (λ > 1) optical bullet is assured for coupling strengths below some critical value. This stability result is also in agreement with the known stability of one dimensional solitons which in this model corresponds to the limit K → 0. So far, we have described the theory on the existence and stability of the LB localized in a few fibers. Now we will show that in fact, the LB can be excited during the propagation from an initial pulse. Specifically, we simulate the propagation of a Gaussian pulse launched in one fiber for the hexagonal geometry: U 0,0 (z = 0, t) = λ k exp(− t 2 8πk 2 λ ), where λ = 14, k = 2.5 (a) and k = 2 (b), so that the energy = 2λ. Figure 9 summarizes the general picture: The left panel highlights the propagation for small λ (small pulse energy), where according to our analysis, stable LBs do not exist. One can observe fast dispersion of the pulse with the energy in the central core vanishing after z = 2, with its energy spreading across the surrounding layers. For large λ (right panel), we observe the formation of LBs localized in the center and the surrounding layer. The observed temporal oscillation can be explained in the following way: The Hamiltonian is conserved in our system and in general the initial value of H and that of the LB are different. The adjustment of the Hamiltonian takes place by emitting radiation in the outside cores which carries the remaining part of the Hamiltonian. Clearly, the strength of these oscillations reflect the efficiency of exciting a LB. For reference, red iso-surfaces are constructed for the central core, and blue ones in the first level cores, represent power values 50 times less that the central power level.
So far we have discussed spatio-temporal light bullet propagation in an optical medium consisting of a multicore fiber or a general array of waveguides which as this and other recent work indicates, can be used in a variety of applications. We have put forward a solid theoretical framework based on asymptotics and stability conditions that not only validates the existence of stable light bullets, representing light localization in space and time, but it demonstrates light bullets are not only achieved for a particular array configuration and instead are generic and stable for a multitude of geometries. We also found that they can be formed as a result of the evolution of a sufficiently intense initial pulse launched into the array. We anticipate our results present an opportunity to design array configurations whose topology is suitable for a particular application including preparing a LB for propagation in bulk Kerr-media [118] or in the atmosphere. Our theoretical approach should explain recent work in nonlinear arrays when the model is extended by incorporating higher order temporal effects [119], or in particular, provide an explanation of the observed additional (large) λ dependent time shift shown in figure 7 of [105]. Clearly such dependence would be quite useful for time delay lines and for efficient coherent pulse combining [120]. The second extension that comes to mind is coupling of active fibers, where at first approximation we can add to the model linear gain and saturation. Finally, while we concentrated our work to the study of light bullets in the anomalous regime, a natural extension is to perform similar studies to discrete spatio-temporal vortices and other nonlinear modes. With respect to the normal dispersion regime, recent experimental results [121] demonstrate the existence of X-waves in 1d-semiconductor waveguide arrays. Theoretically, X-waves carry infinite energy, so once a proper renormalization of the power and the Hamiltonian is done to the approach presented here, an improved existence and stability analysis to 1d and 2d arrays in the normal dispersion regime can be implemented. We conclude this section by presenting detailed results in two configurations [44].

Compression and energy combining in 7-core ring MCF
Consider the ring core geometry. Taking into account the interaction between nearest neighbors only, we can simplify the analysis with the following assumption and neglect all other coupling terms. The dimensionless equations read: The analysis of pulse compression and combining in a MCF with 1D location of the cores was carried out for the initial conditions given by the same Gaussian pulses in each core: where M is the modulation coefficient, and m = 1, ..., N.
The modulation M cos (2πm/N) produces the smooth initial intensity modulation with a maximum in the designated core. This modulation accelerates the compression, and most importantly, makes it less sensitive to the perturbation of the initial parameters. On the other hand, the perturbation must be small enough for the efficient utilisation of each core (M << 1). Based on the results of the previous study [41], we set M = 0.3.
Performing simulations of light pulse propagation initially distributed according to Eq. (28) aimed at the optimisation of pulse compression and energy combining. The dependencies of the basic compression parameters are obtained for the range of parameters P ∈[ 0.31; 1000], τ ∈[ 0.05; 60] for the 7-core MCF. Numerical simulations of Eq. (27) were implemented using the split-step Fourier method with the Padé approximation with scaling and squaring for the matrix exponential (see [48,122]).
The range of values of the initial pulse parameters P and τ were discretized by 250×250 nodes. For each parameter pair (P, τ ) of the grid, the simulation of the dynamics of optical pulses with the form (28) injected into every core of the appropriate MCF was implemented. We tracked the first peak power maximum of the propagating pulse to get the compression or energy conversion at the minimum possible distance along the fiber. In what follows we call such a distance an "optimal" one. Moreover, we took into account only those maxima at which peak power is increased by more than 0.2N times, where N is the number of cores (N > 5).
This approach was applied to cut-off the situations when the peak power maximum is observed without pulse compression. The value 0.2N was chosen empirically. When the energy E is comparable with E cr , the particular initial distribution is not optimal and it can have a few oscillations before the collapse event (here meant to be a sharp energy localisation into a few cores). In this situation, even when strong compression and combining takes place, the position of peaks is sensitive to the initial variations of the power and pulse duration and, for the purpose of this study, we consider the situation as not practical.
The above procedure eliminates this type of dynamics. The maps of combining and compression performance for the 7-core ring MCF are presented in Fig. 10. In Fig. 10a we see coherent effective and coherent beam combining. For the optimal choice of parameters, 83.5% of the initial energy is combined in one core at the distance 10.36 (in the dimensionless variable z). The combined pulse is smooth, and the side satellite peaks have an intensity of 4% of the peak intensity. The total energy in the wings is about 1.9%. Note that the region of maximal combining lies close to the line H = 0. On the other hand, the map of the compression efficiency is quite different. The blue area denotes pairs of parameters P and τ , for which there is stripe narrowing toward high total energies and is confined by the level L D /L NL ≈ 3000. If L D /L NL > 3000, a large nonlinearity destroys a smooth pulse shape before the compression point. In the region of P and τ values for which the compression occurs, the compression factor of up to 307 times can be obtained. The maximal compression is reached at the distance z = 9.68. The evolution of the peak intensities in the different cores in the case of optimal combining (83.5% of all energy combines in the single core) is presented in Fig. 11. We see that the combined pulse has a smooth temporal shape with low-intensity wings. The length of combining is about 10.36 in dimensionless units, due to the slow initial development. In the final stage, as one can see in Fig. 11a, the intensity growth is about exponential. The pulse compression is modest (about 6 times).
The evolution of the intensity in the different cores for the case of maximal pulse compression is presented in Fig. 12. The pulse is compressed over 300 times yet it is still mostly smooth, with noticeable wings. In this case, only about 36% of the energy it contained in the compressed pulse. The interesting feature of this compression process, is a long initial evolution which is consistent with high sensitivity to the initial conditions.

Compression and energy combining in hexagonal MCF
A comprehensive analysis of pulse compression and combining in MCF with 2D hexagonal location of the cores using the Eq. (1) was carried out for initial pulses having the same Gaussian profile in each core, namely We consider here only 7-core (hexagon plus central core) MCFs. The maps of the pulse combining and compression for 7 cores case are presented in Fig. 13. The blue area denotes pairs of parameters P and τ , for which there is no pulse compression, or (b) the initial pulses (29) with this parameters hich break into filaments or are compressed after the distance along the fiber where the first local maximum of peak power is located. Due to the symmetry of the problem, both combining and compression take place in the (0,0) core. We see that the optimal conditions for combining and compression are very different from the results for the ring cores distribution (Fig. 10). The optimal parameters for the compression are very different from similar combining for the ring MCFs. Efficient combining takes place in a much broader range of parameters in the vertical band, insensitive to the pulse duration. This indicates that collapse takes place mainly in the transverse direction. The efficiency of combining and compression is comparable with that of the ring MCFs, but much less sensitive to variation of the initial parameters. The maximum combining efficiency for a 7-core hexagonal MCF is better and equals 91.6% (Fig. 14) at the distance z = 1.78. In contrast to the ring core configurations, a wide region in the parameter space (P, τ ) exists, where the part of the energy in the central core at the compression moment exceeds 70% of the initial energy E. The presence of this region allows us to obtain well-compressed pulses having most of the total energy E, which is of great practical importance. However, in this case a substantial part of the energy goes into the wings, rather than into the central peak of the compressed pulse. Using the 7-core fiber, a maximal pulse compression up to 256 times can be achieved (Fig. 15). In contrast to the combining pattern, at the point of maximal compression a significant fraction of energy is left in the neighbor cores. The peak power increases greatly at the compression point. The best compression occurs in the case of a high power input pulse. Too much nonlinearity (L D /L NL ≈ 4000) destroys the pulse shape before the compression point and confines the maximal pulse compression factor. On the other hand, it is difficult to define the optimal compression case in the presence of high nonlinearity. The pulse compression factor close to its maximal value can be obtained for different pairs of parameters P, τ and the distance to these compression points decreases with growth of the parameter P. The compression distance is sufficiently small (z ≤ 0.1).

OAM modes and the poynting vector singularities in optical fibers
In the remainder of this paper, we would like to discuss a topic which is not directly related to the discrete optical vortices: singularities of the Poynting vector in continuous media(here we examine the case of optical fibers). However, we believe that the analysis is inherently linked to the field of optical vortices which may stimulate new research directions in the study of optical vortices in photonic lattices and MCFs. It is known that optical vortices characterized by the Hilbert factor exp(imφ), describe helical dislocations in the light fields. Here the integer m is a topological charge and the phase distance between consequent coils of the helicoids equals to 2mπ. The phase distance between the neighboring helicoids is 2π. The topological charge of the optical vortex can be defined as: where S is the phase of the scalar light field in question and integration is performed along the contour surrounding the vortex position. The theory and experiments related to optical vortices in free space have a long history (see, e.g. [10,123] and references therein) and generation and propagation of the vortex beams are quite well studied and understood. At the same time, properties of optical vortices in optical fibers, including microstructured fibers, are less explored and this is an emerging area of research in fiber optics. One of the attractive and promising potential applications of vortex modes with orbital angular momentum [124] (OAM modes) in optical fibers is for data transmission in telecommunications [7,9,125]. In [126] it was demonstrated a terabit -scale high -capacity optical communication via OAM multiplexing. In addition to practical applications, the study of optical vortices, as well as the spin and orbital angular momentum of modes in optical fibers, is interesting from the fundamental science which goes beyond the scope of fiber optics [127][128][129][130].
Recently, OAM modes have attracted interest in nonlinear fiber optics because possible nonlinear frequency generation of OAM modes can be exploited in a range of applications from super -resolution microscopy to quantum hyper -entanglement. In [131] the design of a triangular -lattice annular core photonic crystal fiber was proposed for the stable broadband guided -transmission and supercontinuum generation of optical vortex beams in fibers. Using a specially designed optical fiber, the authors in [132] demonstrated experimentally Raman solitons in OAM modes and the first supercontinuum spanning more than an octave (630 nm to 1430 nm) with the entire spectrum in the same polarization as the OAM state. In [133] the Raman gain of different OAM modes for the pump and signal was experimentally characterized. It was demonstrated that Raman gain among OAM modes is independent of the topological charge of the OAM mode. In [134], the first experiments on nonlinear four wave mixing between OAM modes in an optical fiber have been carried out. These experiments show that the large modal space, as well as spin and OAM conservation rules, enable a high diversity of the phase matching conditions.
It is known that vortices and singularities can also be observed in the internal energy flows of light beams [135][136][137]. In this case, singularities of the internal flows occur where the Poynting vector or its transverse component vanishes [14]. In this section, we will briefly discuss the main properties of OAM modes in optical fibers and singularities in the Poynting vector components of the core modes of microstructured fibers.

Phase singularities and OAM modes in optical waveguides (mathematical description)
To the best of our knowledge the morphology of optical dislocations was first considered in [138]. In this work, the authors performed a detailed analysis of phase dislocations in wave trains and in monochromatic optical fields. They demonstrated that the existence of pure and mixed screw dislocations in optical waves. In later works [139], screw dislocations were called optical vortices. Optical vortices have a helicoidal wavefront or a system of |m| helicoids which are located on the same axis. On this axis of the phase dislocation, there is a phase uncertainty and the optical wave amplitude is equal to zero. This means that there is a zero -intensity line at the axis along the direction of wave propagation. In general, dislocations (optical vortices) are typically observed in a wavefront of optical fields with a complex spatial structure (speckled fields) [140]; however, arrays of several screw dislocations can also exist in regular optical waves. The behavior of vortices embedded in nearly monochromatic optical beams [141], i.e. electromagnetic fields in the optical frequency range with a well -defined direction of propagation and with polarization perpendicular to the beam propagation axis, is described by the paraxial propagation equation for the complex scalar electric field E(x, y, z) [142] : where f (E, E * , x, y, z) is a function which can be nonlinear. Since the amplitude of the scalar field must vanish at the singular point and the field phase has to vary continuously around the vortex center the most studied optical vortex field is proportional to r |m| exp(imφ) with the helical phase given by exp(imφ), where r is a radius -vector, φ is an azimthal angle and m is a topological charge. All the above applies to optical fields in the free space. For optical waveguides, most published work related to OAM modes and their potential applications in telecommunications deal with step -index fibers [143], (the step -index fibers with air core [144,145]) or with the graded index fibers [146]. In addition, the generation and propagation of OAM modes was considered for different types of micro -structured optical fibers, such as twisted photonic crystal fibers [147] and hollow core micro -structured fibers [148,149].
In the case of weakly guiding step -index fibers, the scalar wave equation for the transverse component of the electric field of the core modes e t is: where k = 2π/λ, n 1 is a refractive index of the fiber core and β is the propagation constant of the core mode. Each of the OAM modes can be represented by a linear combination of two LP modes of the weakly guiding step -index fiber as follows [143]: where R m (r) = J m (u mn r)/J m (u mn a) if r < a, where a is the fiber core radius and R m (r) = K m (w mn r)/K m (w mn a) if r > a. The parameters are u mn = a ((kn 1 ) 2 − β 2 mn ) and w mn = a β 2 mn − ((kn 2 ) 2 ), correspondingly. Here n 2 is the refractive index of the cladding, C mn is a complex coefficient, m is the order of Bessel function and n is the number of its root. LP modes by themselves cannot carry the orbital angular momentum of light since their angular dependencies are expressed as cos(mφ) and sin(mφ), but their linear combination with the corresponding phase delay of π/2 leads to the appeareance of an OAM mode with a helical phase given by exp(imφ). OAM modes have a z -projection of the linear density of the orbital angular momentum which is proportional to the index m (topological charge) [150].
In the case when the fiber's refractive index gradient cannot be neglected, solving the full vector equation for the transverse component of the electric field e t is required: OAM modes can also be represented as a π/2 -phase -shifted linear combination of vector modes which are solutions to Eq. (34) [6]. Let us consider the formation of OAM modes for the step -index fiber using the example of a linear combination of HE m+1,n vector modes. The expressions for even and odd HE m+1,n modes in Cartesian coordinates are: HE odd m+1,n = R m (r)( xsin(mφ) + ycos(mφ))exp(iβ mn z) Then, choosing a new coordinate system in the form of σ ± = x ± iy, where σ ± represents circular polarization at which each photon carries spin angular momentum with value of S = ±1, the new mode can be expressed as: It can be seen from (36) that the new state mn (r, φ, z) is an OAM mode with a helical phase given by exp(imφ) and the spin of the OAM modes is aligned with the sign of the toplogical charge of m. OAM modes can be obtained from even and odd EH m−1,n modes in the same way, but with the spin being anti -aligned with the sign of topological charge of m. The non degenerate modes TE 0,n and TM 0,n can also be represented in the same basis [6]. On the experimental side, the authors in [151], present a complete generation and detection scheme for radially modulated vector -vortex modes for step -index fibers, by employing both dynamic and geometric phase control of light.

Features of angular momentum of light and the poynting vector singularities in optical waveguides
General features that include spin and orbital angular momentum were studied in [152] for structured monochromatic optical fields propagating in dispersive inhomogeneous isotropic media. This task is associated with a more general problem known as the Abraham -Minkowski dilemma where one considers kinetic (Poynting -like) versus canonical (spin -orbital) pictures. For this, the authors introduced a novel Minkowski -type momentum, spin and orbital angular momentum densities of the field and it was shown that the canonical momentum density can be associated with the local wavevector (which can be defined as a phase gradient) while the Poynting -Abraham momentum density can be associated with the energy flux and the group velocity of the propagating wave.
In [127], the authors applied these results to describe the angular momentum of light in dielectric step -index fibers. They showed that the total canonical (orbital + spin) angular momentum of the modes in a dielectric step  [153]. The authors demonstrated the formation of singularities in the Poynting vector transverse component for a fiber with a parabolic refractive index profile. In this case, the axial projection of the angular momentum of the vortex mode was determined according to an approach related to the kinetic Abraham momentum density. The singularities of the Poynting vector transverse component of the core modes were determined using the following equations P x (x, y) = P y (x, y) = 0, where P x , P y are the projections of the Poynting vector transverse component on the Cartesian coordinate axis and (x, y) are the coordinates of the same point in the fiber cross-section.
Based on a hydrodynamic approach [154], the pattern of stream lines of the Poynting vector transverse component of the leaky fundamental mode of hollow core microstructured optical fiber was used to study the properties of radiation outflow from the air core. It was shown that the stream lines of the transverse Poynting vector component of the fundamental air -core modes in a capillary, polygonal hollow core waveguides and in the negative curvature hollow core fibers [155], form vortices with the centers located on the fiber axis and in the cladding elements. This process is schematically shown in Fig. 16 for a capillary and for the negative curvature hollow core fiber. Based on the obtained results, it was assumed that vortices arising in the Poynting vector transverse component of the core modes should affect the loss level in hollow core microstructured fibers.
In [156], the authors showed that vortices in the Poynting vector transverse component of the fundamental mode can form both in the hollow core micro -structured fibers and in all solid band gap fibers [157]. A study of the loss level of the fundamental core modes and the structure, shape, and arrangement of singularities of their Poynting vector transverse component suggest a close relationship between them. This can be well demonstrated in the one layer all solid band gap fiber (Fig. 17). The geometric and optical parameters of the fiber are as follows: pitch = 12 μm, ratio d/ = 0.33 (d is the cladding rod diameter), the refractive index of the surrounding glass matrix n g = 1.45 and the refractive index of the cladding rods n r = 1.5. Large losses appear in the transmission band centered at a wavelength of 1.5 μm. Losses were approximately two orders of magnitude lower in the transmission band centered at a wavelength of 1 μm. As can be seen from Figs. 17, 18, and 19, the structure of the radial component of the Poynting vector of the fundamental core mode, as well as the structure of stream lines and vortex centers, are completely different for the cases of small and large losses. This suggests that the vortex structure formed in the fiber cladding determines the trajectories along which the mode energy flows out of the core.
The problem associated with the study of singularities and vortices arising in the Poynting vector of the core modes of micro -structured fibers is at the very initial stages and it seems clear that it can lead in the future to interesting results from both fundamental (the Abraham -Minkowski dilemma) and practical points of view (reduction of loss level by orders of magnitude in microstructured fibers).

Conclusions
In this review we discussed some of recent advances in the field of discrete vortex systems and numerical methods for modeling of light propagation in multi-core fibers. This area is growing fast and we did not intend here to give a comprehensive overview of all interesting works in the field; for example we did not address studies of quantum vortices. Instead, our focus here was on topics close to our research interests. With this in mind, we would like to summarise the key points discussed above.
• Light in discrete or discrete-continuous systems can form rotating structure -discrete vortices having unique topological properties. • There is no general exact analytical solutions describing nonlinear discrete vortices. However, in some situations, approximate solutions can describe the local structure of the discrete vortices. For instance, in the case of MCF with a small number of cores and assuming rotational symmetry, one can use two-mode approximation to describe such vortex structures. • Discrete vortices are typically more robust against azimuthal perturbations than continuous vortex structues. • Stable discrete vortices in interacting nonlinear multicore fiber systems with one central and several periphery cores have been shown to exist. Although the presence of the central core and nonlinearity generally tends to destabilize vortices, the coherent propagation of steady nonlinear vortices with high different regimes of the core mode energy leakage and accordingly, to different loss levels. • The study of vortices properties arising in the Poynting vector transverse component of the core modes can make it possible to find the optimal geometric structures of the micro -structured fibers from the point of view of achieving a minimum loss for the core modes for given material parameters of the fiber.
Control of the topological charge of optical vortices can open a new area of applications [158]. Optical vortices are naturally linked to topological photonics and we anticipate that this field will be further advancing both on fundamentals and in applications. Fig. 19 The streamlines of the transverse component of the Poynting vector (grey lines) and the lines of zero intensity P t = (P x = 0(redlines); P y = 0(greenlines)) for all solid band gap fiber with a cross section shown in Fig. 17 at a wavelength of λ = 1.5μm.
(2021) 17:23 Page 24 of 28 each fixed layer n on an evolutional variable z. Specifically, we define the sequence V 0 , V 1 , · · · , V s (s is the iteration number), that converges to the solution U n+1 on the (n + 1)-th layer. Finally, after an initial approximation for V k,0 , k = 1, ..., N, which can be found, for example, by explicit scheme or from the previous step by the spatial variable z, one can iterate the linearized equation [47] iI − β 2 αh 2 to obtain the series V k,0 , V k,1 , · · · , V k,s , until |V k,s − V k,s−1 | < ε for every k = 1, ..., N, where I is the identity operator and ε ∼ 10 −8 is the tolerance level. The upper limit of iterations is equaled to 30. The Eq. (43) is obtained by letting U n+1 k = V k,s+1 in Eq. (41). The dissipative parameter c was chosen empirically and was set to 0.05 to minimize the error of calculations. We implemented the periodic boundary conditions by using the Sherman-Morrison formula [160,161].

Split-step fourier method with the padé approximation of the matrix exponential
This numerical algorithm is a spectral numerical method with high accuracy in the calculation by the time variable [46,48]. The split-step Fourier method has been adapted to a system of coupled NLSE in a special way [35,162]. We first rewrite the system (38) in the operator form: where U = (U 1 , ..., U N ) T is a vector of complex variables, D is a linear differential operator defining dispersion and linear couplings among variables U k , k = 1...N, andN is a nonlinear operator. These operators have the following structure for the system of equations (38) considered above: where I is an identity operator, and C is a matrix of linear couplings. An approximate solution in SSFM is obtained, assuming that during the propagation of the optical field on one step h of the evolutional variable z, linear and nonlinear effects act independently. In particular, the propagation from z to z + h is carried out in two steps. In the first step, only nonlinearity acts and the operatorD = 0. In the second step, only linear effects are taken into account, and therefore the nonlinear operatorN = 0. In practice, the following symmetrized form of the algorithm is used to achieve the second order approximation: When SSFM is generalized on the NLSE system (44), the calculation of the nonlinear diagonal operator exp(αhN), α = {0.5, −0.5} is split by calculation of the individual actions of the nonlinear operator for each variable U k , k = 1, ..., N, using an explicit algorithm [163]. The linear oper-atorD is a non diagonal matrix, therefore, calculation of the matrix exponential exp(hD) is needed for implementation of linear step. The following scheme was used to calculate it. The matrix exponential exp(hD) is calculated in Fourier space, using the representation exp(hD)Ũ(z, t) = F −1 t exp[ hD(−iω)] F tŨ (z, t), (49) where where z = nh, and the operatorD(−iω) was obtained from the dispersion operator (45) by replacing derivative ∂/∂t by −iω, where ω is the frequency in Fourier space. The resulting matrix exp[ hD(−iω)] is constant, so in the case of a uniform grid in z it is sufficient to carry out its calculation once at the beginning in the most accurate way. In practice, the numerical Padé approximation is the most used method for computing a matrix exponential [122,164] . The Padé approximant for a matrix exponential can be written in the form