Optical vortices in waveguides with discrete and continuous rotational symmetry

Abstract

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–10]). Vorticity is characterized by an integer winding number (topological charge), which is 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–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–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–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–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–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–48] or the nonlinear complex Ginzburg-Landau equation [49–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:

$$ i\frac{\partial}{\partial z}U_{nm} + (\underline{\underline{C}} \ \underline{U})_{n,m} + \frac{\partial^{2}}{{\partial t^{2}}}U_{nm} + 2|U_{nm}|^{2}U_{nm} =0, $$

(1)

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 - $(\underline {\underline {C}}\ \underline {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:

$$\begin{array}{@{}rcl@{}} (\underline{\underline{C}}\ \underline{U})_{nm}^{square} =c(U_{n-1,m} + U_{n+1,m} + U_{n,m-1} + U_{n,m+1}) \\ (\underline{\underline{C}}\ \underline{U})_{nm}^{hexagon}\! = \!c(U_{n-1,m} \!+ \!U_{n+1,m} \!+ \!U_{n,m-1} \!+ \!U_{n,m+1} \!+ \\ U_{n,m+1} + U_{n+1,m+1} + U_{n+1,m-1}) \end{array} $$

(2)

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 waveguides. Finally, an “Appendix describes the computational tools used in the study of vortex dynamics in MCFs.

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–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–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:

$$\begin{array}{*{20}l} i\frac{\partial U_{n,m}}{\partial z} +c (U_{n-1,m} + U_{n+1,m} + U_{n,m-1} \\ + U_{n,m+1} - 4U_{n,m}) + |U_{n,m}|^{2}U_{n,m} = 0 \end{array} $$

(3)

derived from the corresponding Hamiltonian

$$ {}H \!= \!\sum_{n,m} c |U_{n+1,m} - U_{n},m|^{2} +c |U_{n,m+1} - U_{n},m|^{2} - \frac{1}{2}|U_{n,m}|^{4} $$

(4)

The system also conserves total power (also called in the literature the particle number) $N = \sum _{n,m} |U_{n,m}|^{2}$. Note that Eq. 3 has a family of stable two-dimensional discrete solitons (see, e.g. [64–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^{(0)}_{n,m}= Ae^{i(z+ \theta _{nm}-\theta _{0})}$, where C satisfies (1−A²)A=0, with two possible solutions A=0,A=1 or

$$ U^{(0)}_{n,m} = \left\{ \begin{array}{cc} e^{i(z+ \theta_{nm}-\theta_{0})} & (n,m) \in S\\ 0 & (n,m) \not\in S \end{array}\right. $$

(5)

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, θ₁,θ₂,....,θ_N+1=θ₁, the condition for persistence for c>0 small, based on the IFT is the “zero power flow”

$$ \sin(\theta_{j} - \theta_{j+1}) + \sin(\theta_{j} - \theta_{j-1}) = 0, \ \ 1 \le j \le N $$

(6)

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 $\theta _{1} = 0, \ \theta _{2} = \frac {\pi }{2}, \ \theta _{3} = \pi, \ \theta _{4} = \frac {3\pi }{2}$. A second example for 8 sites is a one parameter family:

$$ \begin{array} {cc} \theta_{1} = 0, \theta_{2} = \theta, \theta_{3} = 2\theta, \theta_{4} = 3\theta, \ \ \theta_{5} = 4\theta, \theta_{6} = \pi + \\ 3\theta, \theta_{7} = 2\pi + 2\theta, \theta_{8} = 3\pi + \theta \end{array} $$

(7)

with charge = 1 for $\theta = \frac {\pi }{4}$, charge = 2 for $\theta = \frac {\pi }{2}$ and charge = 3 for $\theta = \frac {3\pi }{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:

$$ i\frac{\partial U_{n}}{\partial z} + c (e^{i\phi}U_{n-1} + e^{-i\phi}U_{n+1}) + \gamma |U_{n}|^{2}U_{n} = 0 $$

(8)

U_n+N=U_n. Here the phase ϕ is proportional to the twist rate. To illustrate the topological nature of the twist, consider the N=4 array with $\phi =\frac {\pi }{4}$ in the linear case γ=0. An exact solution reads

$$\begin{array}{@{}rcl@{}} U_{1} &=& \cos(\sqrt{2}c z),\\ U_{2} &=& \frac{1}{2}(i-1)\sin(\sqrt{2}c z), \\ U_{3} &=& 0,\\ U_{4} &=& \frac{1}{2}(i+1)\sin(\sqrt{2}c z), \end{array} $$

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 $\phi = \frac {\pi }{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–83] with a refractive index profile along the radial direction in a form of zero order or higher order Bessel functions $R(\eta,\xi,) =J_{n}^{2} (2b_{lin})^{1/2} (\eta ^{2}+\xi ^{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–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–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 [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–88], coherent beam combining and new concepts and techniques in the field of high-power fibre lasers [91–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–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 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 low-dimensional version of the discrete nonlinear complex Ginzburg-Landau equation in a weak coupling approximation, assuming significant interactions between nearest neighboring cores [40, 84]:

$$\begin{array}{@{}rcl@{}} i\frac{d U_{0}}{dz}+\beta_{0}U_{0}+\sum_{m=1}^{M}C_{0,m} U_{m} +2\gamma_{0}\left|U_{0}\right|^{2}U_{0} \\ = i\left(G\left(U_{0}\right)-\alpha_{0}\right)\frac{U_{0}}{2}, i \frac{{dU}_{m}}{dz}+\beta_{m}U_{m} + \\ C_{0,m}U_{0}+C_{m,m+1}U_{m+1}+C_{m,m-1}U_{m-1} + \\ 2\gamma_{m}\left|U_{m}\right|^{2}U_{m} =i\left(G(U_{m})-\alpha_{1}\right)\frac{U_{m}}{2}. \end{array} $$

(9)

Fig. 1

Schematic of an MCF with M=5 periphery cores and the central core. Geometric parameters that dominantly determine coupling coefficients C₀ and C₁, 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

Here, U₀ 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 β₀ 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 γ₀ 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=β₁, C_m,m−1=C_m,m+1=C₁, and γ_m=γ₁ and identical coupling between central core and periphery cores C_0,m=C₀ the Eq. 9 can be rewritten as:

$$\begin{array}{@{}rcl@{}} i\frac{d U_{0}}{dz}+\beta_{0}U_{0}+C_{0}\sum_{m=1}^{M}U_{m} + \\ 2\gamma_{0}\left|U_{0}\right|^{2}U_{0} = 0, i \frac{{dU}_{m}}{dz} + \\ \beta_{1}U_{m}+C_{0}U_{0}+C_{1} (U_{m+1} + \\ U_{m-1}) +2\gamma_{1}\left|U_{m}\right|^{2}U_{m}=0. \end{array} $$

(10)

Two conserved quantities characterize light propagation through the MCF system (10): the norm (total power) P and the Hamiltonian H:

$$\begin{array}{*{20}l} P&=&\left|U_{0}\right|^{2}+\sum_{m=1}^{M}\left|U_{m}\right|^{2}, \end{array} $$

(11)

$$\begin{array}{*{20}l} H&=&-\sum_{m=1}^{M}\left[-\beta_{1}\left|U_{m}\right|^{2}-C_{1} \left(U_{m+1}U_{m}^{*}+U_{m+1}^{*}U_{m}\right) \right. \\ &-& \left. C_{0}\left(U_{0}U_{m}^{*}+U_{0}^{*}U_{m}\right) -\gamma_{1}\left|U_{m}\right|^{4} \right] \\ &-& \beta_{0}\left|U_{0}\right|^{2}-\gamma_{0}\left|U_{0}\right|^{4}, \end{array} $$

(12)

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:

$$ \frac{d \delta{U}}{dz}=i\hat{M}\delta{U} $$

(13)

where δU is a column matrix of small perturbations $\delta {U}=[\delta {U}_{0}, \delta {U}_{1},... \delta {U}_{M}, \delta {U}_{0}^{*}, \delta {U}_{1}^{*}, \delta {U}_{M}^{*}],$ and $\hat {M}$ is the corresponding eigenvalue matrix with dimensions 2N×2N 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 matrix $\hat {M}$ 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):

$$ \lambda_{k} =\beta_{1} C_{1} \cos{(2k\pi/M)} $$

(14)

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

$$\begin{array}{*{20}l} \lambda_{M,M+1} =\beta_{1R}+ C_{1R} +i \frac{\beta_{1I}+2 C_{1I}+ \beta_{0I}}{2} \\ \pm \sqrt{{MC}_{0}+\left(C_{1R}+ \frac{\beta_{1I}+2 C_{1I} - \beta_{0I}}{2}\right)} \end{array} $$

(15)

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)

$$\begin{array}{@{}rcl@{}} u_{n}(z)&=&a_{M}u^{(M)}e^{i\lambda_{M} z}+a_{M+1}u^{(M+1)}e^{i\lambda_{M+1}z} \\ &+& \frac{1}{\sqrt{M}}\sum_{k=1}^{M-1}a_{k}e^{i 2\pi k n /M +i\lambda_{k} z}\\ u_{0}(z)&=&b_{M}U^{(M)}e^{i\lambda_{M} z}\\ &+& b_{M+1}U^{(M+1)}e^{i\lambda_{M+1} z},\, n=1,...,M, \end{array} $$

(16)

where

$$\begin{array}{@{}rcl@{}} u^{(j)}&=&\frac{\lambda_{j}+i\beta_{0I}}{\sqrt{M^{2}C_{0}^{2}+M|\lambda_{j}+i\beta_{0I}|^{2}}},\\ U^{(j)}&=&\sqrt{\frac{{MC}_{0}^{2}}{{MC}_{0}^{2}+|\lambda_{j}+i\beta_{0I}|^{2}}},\quad j=M,M+1 \end{array} $$

(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

$$ {}U_{0}=0,\qquad U_{m}=A\, e^{i\frac{2\pi S}{M}m}=A\, e^{i\alpha m},\quad m=1,...,M. $$

(18)

Note that $\sum _{m=1}^{M} \exp [i 2\pi S m/M]=0$ for any integer S. The vortex solutions are eigenmodes with arbitrary amplitude associated with the eigenvalues

$$ {}\lambda_{k}=\beta_{1R}+2C_{1R} \cos{\left(\alpha+\frac{2\pi k}{M}\right)},\, k = 1,2,...,M-S. $$

(19)

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

$$\begin{array}{*{20}l} A^{2} \ge \frac{C_{1}}{2\gamma_{1}} \cos(\alpha)[1-\cos(2\pi j S/M)] \quad \\ {\text{and}} \quad \frac{4p+1}{4} \le \frac{S}{M} \le \frac{4p+3}{4} \end{array} $$

(20)

where index j denotes different eigenvalues (j=0,...,M−1), and p is an integer. The following sets of stable vortices with unlimited power are possible: S=1, M=3,4; S=2,M=3,4,5,6,7,8; S=3, M=4,5,6,7,8,9,10,11,12; S=4, M=3,6,7,…,16; and so on.

II) If cos(α)≥0 then stable vortices are characterized by

$$\begin{array}{*{20}l} A^{2} \le \frac{C_{1}}{2\gamma_{1}} \cos(\alpha)[1-\cos(2\pi j S/M)] \quad \\ {\text{and}} \quad \frac{4p-1}{4} \le \frac{S}{M} \le \frac{4p+1}{4} \end{array} $$

(21)

This inequality gives the following sets of stable vortices in the low power region: S=1, M=5,6,7,...; S=2, M=9,10,...; S=3, M=3,13,14,...; S=4, M=4,5,17,18,...; and so on.

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.

Fig. 2

Schematic of the stability regions for vortices with topological charge S=1 (a), S=2 (b), S=3 (c), and S=4 (d) with respect to the number of cores in the MCF system. Blue and green bars in plots denote the parameter areas with stable vortices in the MCF without and with central core, respectively

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.

Fig. 3

The examples of the vortex with S=2 propagation in the MCF with M=6 periphery cores. The phase evolution of dynamically stable vortex with P=60.03 (left), and dynamically unstable vortex with P=10.8 (right)

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–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–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\frac {\partial }{\partial z} U(x,y, z) + \Delta _{T} U+|U|^{2\sigma }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) that are transverse to the direction of propagation z. At the so called critical case (σd=2, d=transverse dimension), collapse is easily arrested by various additional terms in the model even if they are small. Most common examples are nonlinear losses and nonlinearity saturation.

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–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)=\sum _{nm}U_{nm}(z,t)F(x-x_{m},y-y_{m})e^{i(kz-\omega 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:

$$ i\frac{\partial}{\partial z}U_{nm} + (\underline{\underline{C}} \ \underline{U})_{n,m} + \frac{\partial^{2}}{{\partial t^{2}}}U_{nm} + 2|U_{nm}|^{2}U_{nm} =0, $$

(22)

where $(\underline {\underline {C}}\ \underline {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 coupling operators are:

$$\begin{array}{@{}rcl@{}} (\underline{\underline{C}}\ \underline{U})_{nm}^{square} =c\left(U_{n-1,m} + U_{n+1,m} + U_{n,m-1} + \right. \\ \left. U_{n,m+1}\right) \\ (\underline{\underline{C}}\ \underline{U})_{nm}^{hexagon} =c\left(U_{n-1,m} + U_{n+1,m} + U_{n,m-1} + U_{n,m+1} + \right. \\ \left. U_{n,m+1} + U_{n+1,m+1} + U_{n+1,m-1}\right) \end{array} $$

Fig. 4

Square and hexagonal waveguide structures

This system has two known conserved quantities; the Hamiltonian

$H=\sum _{nm} \int (N(\underline {U}; \underline {U})_{nm} + |\frac {\partial U_{nm}}{\partial t}|^{2}- |U_{nm}|^{4})dt$ and the total power $P=\sum _{nm}\int |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:

$$\begin{array}{@{}rcl@{}} -\frac{\delta(H+\lambda P)}{\delta A^{*}_{nm}}= -\lambda A_{nm} + (\underline{\underline{C}} \ \underline{A})_{nm} \\ + \frac{\partial^{2}}{{\partial t^{2}}}A_{nm} + 2|A_{nm}|^{2}A_{nm} = 0 \ \end{array} $$

(23)

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_{nm}^{(0)}(t;\lambda)+ A_{nm}^{(1)}(t;\lambda) + A_{nm}^{(2)}(t;\lambda) +..$. The analytical derivation of the asymptotic solutions is straightforward and we only present here the expressions for the central core and the first layer: The 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

$$\begin{array}{@{}rcl@{}} A_{0,0}(t)=a_{0}(t) + O(\frac{1}{\lambda^{3/2}}) = \frac{\sqrt{\lambda}}{cosh(\sqrt{\lambda} t)}+O\left(\frac{1}{\lambda^{3/2}}\right) \\ A_{\pm 1,0}^{square}(t) = A_{0,\pm 1}^{square} = a_{1}(t) +\\ O\left(\frac{1}{\lambda^{5/2}}\right) = \frac{c}{2\sqrt{\lambda}} [ e^{\sqrt{\lambda} t} ln(1+e^{-2\sqrt{\lambda} t}) +\\ e^{-\sqrt{\lambda} t} ln(1+e^{2\sqrt{\lambda} t}) ] + O\left(\frac{1}{\lambda^{5/2}}\right) \\ A_{0,\pm 2}^{hex}(t) = A_{1,\pm 1}^{hex} = a_{1}(t) + O\left(\frac{1}{\lambda^{3/2}}\right) \\ = \frac{c}{2\sqrt{\lambda}} [ e^{\sqrt{\lambda} t} ln(1+e^{-2\sqrt{\lambda} t}) \\ + e^{-\sqrt{\lambda} t} ln(1+e^{2\sqrt{\lambda} t}) ] + O\left(\frac{1}{\lambda^{3/2}}\right) \end{array} $$

It should be stated that these solutions fail to be uniformly valid beyond $|t| \ge O(\sqrt {\lambda })$. In fact at the pulse tails, all waveguides have solutions of the same order $A_{nm}=O(t e^{-\sqrt {\lambda } |t|}), t \ge \sqrt {\lambda }$. 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₀(t),a₁(t) coincide with numerical solutions A₀₀(t) and A_±1,0(t)=A_0,±1(t) of Eq. (23) up to the order $O(\frac {1}{\lambda ^{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.

Fig. 5

Peak power (numerical (color) and theory (black)) vs the parameter λ in central (a) and closer neighbor waveguides (b) for 33x33 square and 65x33 hexagonal waveguides structures

Fig. 6

Comparison of numerical solutions (color) with their analytical approximations (black) for different λ in the central waveguide (top) and closer neighbor waveguide (bottom) for 33x33 square (green) and 65x33 hexagonal (red) waveguides structures

Fig. 7

3D power and phase time distributions vs parameter λ for the central (left) and closer neighbor waveguide (right) for a square waveguides structure

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].

Fig. 8

From left to right: (a) dependence of the total power on the parameter λ; (b) dependence of the Hamiltonian on the parameter λ; (c) the Hamiltonian vs the total power for 33x33 square (green) and 65x33 hexagonal (red) waveguides structures

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

$$\begin{array}{*{20}l} -\frac{\partial g_{nm}}{\partial z}=-(\underline{\underline C}\underline{f})_{nm} - \frac{\partial^{2} f_{nm}}{\partial {t^{2}}} \\ +\lambda f_{nm}-6|A_{nm}|^{2}f_{nm} = (H_{-}f)_{nm} \end{array} $$

(24)

$$\begin{array}{*{20}l} \frac{\partial f_{nm}}{\partial z}=-(\underline{\underline C}\underline{g})_{nm} - \frac{\partial^{2} g_{nm}}{\partial {t^{2}}} \\ +\lambda f_{nm}-2|A_{nm}|^{2}g_{nm} = (H_{+}g)_{nm} \end{array} $$

(25)

or $-\frac {\partial ^{2} f_{nm}}{\partial {z^{2}}}= [ H_{+}(H_{-}f) ]_{nm}$. Defining inner products for square arrays:

$$\begin{array}{@{}rcl@{}} < \underline{f},H_{+}\underline{f} > = \sum_{nm}\int f_{nm} \left(H_{+}f_{nm} \right)dt = \\ \sum_{nm}\int [|A_{nm}|^{2} \left(\frac{\partial }{\partial t}\frac{f_{nm}} {A_{nm}}\right)^{2} ] dt + \sum_{nm}\int [ (\sqrt{\frac{A_{n,m-1}}{A_{nm}}}f_{nm} - \\ \sqrt{\frac{A_{n,m}}{A_{n,m-1}}}f_{n,m-1})^{2} ] dt + \sum_{nm}\int [ (\sqrt{\frac{A_{n-1,m}}{A_{nm}}}f_{nm} - \\ \sqrt{\frac{A_{n,m}}{A_{n-1,m}}}f_{n-1,m})^{2} ]dt \end{array} $$

and for hexagonal arrays:

$$\begin{array}{@{}rcl@{}} <\underline{f},H_{+}\underline{f} >= \sum_{nm}\int f_{nm}(H_{+}f_{nm})dt = \\ \sum_{nm}\int [ |A_{nm}|^{2} \left(\frac{\partial }{\partial t}\frac{f_{nm}} {A_{nm}}\right)^{2} ]dt + \sum_{nm}\int [ (\sqrt{\frac{A_{n,m-1}}{A_{nm}}}f_{nm} - \\ \sqrt{\frac{A_{n,m}}{A_{n,m-1}}}f_{n,m-1})^{2} ]dt + \sum_{nm}\int [ (\sqrt{\frac{A_{n-1,m}}{A_{nm}}}f_{nm} - \\ \sqrt{\frac{A_{n,m}}{A_{n-1,m}}}f_{n-1,m})^{2} ]dt + \sum_{nm}\int [ (\sqrt{\frac{A_{n+1,m+1}}{A_{nm}}}f_{nm} - \\ \sqrt{\frac{A_{n,m}}{A_{n+1,m+1}}}f_{n+1,m+1})^{2} ] dt + \sum_{nm}\int [(\sqrt{\frac{A_{n+1,m-1}}{A_{nm}}}f_{nm} - \\ \sqrt{\frac{A_{n,m}}{A_{n+1,m-1}}}f_{n+1,m-1})^{2} ] dt \end{array} $$

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 $\underline {F}$ for which $< -\underline {F},H_{+}\underline {F} > \ < 0, \ < \underline {F}, \underline {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 $\frac {d}{d\lambda }P= \frac {d}{d\lambda } < \underline {A},\underline {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 $\lambda _{c}=(\frac {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)= \sqrt {\frac {\lambda }{k}}exp(-\frac {t^{2}}{8\pi k^{2}\lambda })$, 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.

Fig. 9

Power iso-surfaces of weak linear (a) and powerful nonlinear (b) input Gaussian pulses with the same energy during propagating along the hexagonal waveguide geometry

So far we have discussed spatio-temporal light bullet propagation in an optical medium consisting of a multi-core 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

$$ C_{m,m+1} = C>0 $$

(26)

and neglect all other coupling terms. The dimensionless equations read:

$$ i\frac{\partial U_{m}}{\partial z} = -\frac{\partial^{2} U_{m}}{\partial t^{2}} - (U_{m+1} - 2U_{m} + U_{m-1}) - |U_{m}|^{2} U_{m}. $$

(27)

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:

$$ U_{m}(0,t) = \sqrt{P} \exp\left(-\frac{t^{2}}{2 \tau^{2}}\right) \left[ 1 + M\cos\left(\frac{2 \pi m}{N}\right) \right], $$

(28)

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 no pulse compression; instead the initial pulses (28) break into clusters as a result of the modulation instability. It is interesting to note that isolines of pulse compression factor in Fig. 10b correlate with isolines of the ratio of the dispersion length L_D=τ²/|β₂| and the nonlinear length L_NL=1/(γP). The zone of optimal compression is the 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.

Fig. 10

Dependence of the percentage of total energy E combined in the core with pulse compression (a) and the pulse width compression factor in the logarithmic scale (b) on the initial pulses [Eq. (28)] parameters P and τ for 7-core MCF with ring core configuration and modulation coefficient M=0.3. The isoline E=E_cr=4π is indicated by the yellow dots, the level H=0 is depicted by white dashes

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).

Fig. 11

Evolution of the input Gaussian pulses [Eq. (28)] with the parameters P=0.5, τ=2.0 and M=0.3 injected in all cores of ring 7-core MCF (best combining) (a). Corresponding input intensity distribution in core 7 (dashed red line) and the distribution at compression point (solid blue line) are shown in logarithmic (b) and normal (c) scales. 83.5% of total input energy E is combined in the core 7. The pulse width (FWHM) is reduced by 5.74x. The peak power increases in 19.5x

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.

Fig. 12

Evolution of the input Gaussian pulses [Eq. (28)] with the parameters P=0.436, τ=42.5 and M=0.3 injected in all cores of ring 7-core MCF (best compression) (a). Corresponding input intensity distribution in core 7 (dashed red line) and the distribution at compression point (solid blue line) are shown in logarithmic (b) and normal (c) scales. 35.8% of total input energy E is combined in the 7th core. The pulse width (FWHM) is reduced by 307.6x. The peak power increases in 123.7x

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

$$ U_{n,m}(z=0,t) = \sqrt{P} \exp\left(-\frac{t^{2}}{2 \tau^{2}}\right). $$

(29)

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).

Fig. 13

Dependence of the percent of total energy E combined in the core with pulse compression (a) and the pulse width compression factor in the logarithmic scale (b) on the initial pulses [Eq. (29)] parameters P and τ for 7-core hexagonal MCF without modulation. The level H=H_c is depicted by white dashes

Fig. 14

Evolution of the input Gaussian pulses [Eq. (29)] with the parameters P=0.687 and τ=1.775 injected in all cores of hexagonal 7-core MCF (best combining) (a). The input intensity distribution in the central core (dashed red line) and the distribution at compression point (solid blue line) are shown in logarithmic (b) and normal (c) scales. 91.6% of total input energy E is combined in the core (0,0). The pulse width (FWHM) is reduced by 6.37x. The peak power increases in 37.7x

Fig. 15

Evolution of the input Gaussian pulses [Eq. (29)] with the parameters P=536 and τ=3.2 injected in all cores of a hexagonal 7-core MCF (best compression) in the logarithmic scale (a). The input intensity distribution in the central core (dashed red line) and the distribution at compression point (solid blue line) are shown in logarithmic (b) and normal (c) scales. 16.6% of total input energy E is combined in the core (0,0). The pulse width (FWHM) is reduced by 256.3x. The peak power increases in 21.85x

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:

$$ \\ m = \frac{1}{2\pi}\oint{\vec{{\nabla}}S{d\vec{l}}} $$

(30)

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–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–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] :

$$ i\frac{{\partial}E}{{\partial}z} +\frac{1}{2}{{\nabla}^{2}E}+f(E,E^{*},x,y,z) = 0 $$

(31)

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:

$$ {{\nabla_{t}^{2}}{e_{t}}}+\left(k^{2}n_{1}^{2}-\beta^{2}\right){e_{t}} = 0\\ $$

(32)

where k=2π/λ, n₁ 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]:

$$ \Psi_{mn} (r, \phi, z) = C_{mn} R_{m}(r) exp(im\phi) exp(i\beta_{mn} z)\\ $$

(33)

where R_m(r)=J_m(u_mnr)/J_m(u_mna) if r<a, where a is the fiber core radius and R_m(r)=K_m(w_mnr)/K_m(w_mna) if r>a. The parameters are $u_{mn} = a \sqrt {(({kn}_{1})^{2}-\beta _{mn}^{2})}$ and $w_{mn} = a\sqrt {\beta _{mn}^{2}-(({kn}_{2})^{2})}$, correspondingly. Here n₂ 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 ${\vec {e_{t}}}$ is required:

$$ {{\nabla^{2}}\vec{e_{t}}}+(k^{2}n_{1}^{2}-\beta^{2})\vec{e_{t}} +{\nabla_{t}}(\vec{e_{t}}*ln{\nabla}n_{1}^{2})= 0\\ $$

(34)

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}_{m+1,n}^{even} = R_{m}(r)({\vec{x}cos(m\phi)-\vec{y}sin(m\phi)})exp(i\beta_{mn}z) $$

(35)

$$ {HE}_{m+1,n}^{odd} = R_{m}(r)({\vec{x}sin(m\phi)+\vec{y}cos(m\phi)})exp(i\beta_{mn}z)\\ $$

(36)

Then, choosing a new coordinate system in the form of $\vec \sigma ^{\pm } = {\vec {x}\pm \vec {iy}}$, where $\vec \sigma ^{\pm }$ represents circular polarization at which each photon carries spin angular momentum with value of S=±1, the new mode can be expressed as:

$$ {}\vec{\Psi}_{mn}(r,\phi,z) = {HE}_{m+1,n}^{even}\pm {iHE}_{m+1,n}^{odd}= R_{m}(r)\vec{\sigma}^{\pm} e^{\pm{i}m\phi}e^{i\beta_{mn}z}\\ $$

(37)

It can be seen from (36) that the new state $\vec {\Psi }_{mn}(r,\phi,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 - index fiber always equals integer values m. The integer value m is the topological charge of the longitudinal mode field component in units of h (Planck’s constant) per photon. In this case, the transverse components of the mode fields are given in terms of the basis of circular polarization attached to the transverse Cartesian coordinates (see “Phase singularities and OAM modes in optical waveguides (mathematical description)” subsection). Also, it was shown that the helicity of guided modes for the dielectric step - index fiber characterizes the degree of chirality of the electromagnetic fields. The helicity of the guided modes in contrast to spin angular momentum, can take any values in the range of [-1; 1] in units of h per photon.

The first study of stream lines of the Poynting vector transverse component of fiber core modes was carried out in [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 micro - structured 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 micro - structured fibers.

Fig. 16

Vortex – supported guidance in micro – structured fibers. (left) Monochromatic radiation in the leaky core mode of the hollow capillary. The axial component of the Poynting vector is from left to right (wavy red arrows). The mode losses are determined by the transverse component of the Poynting vector whose streamlines are represented by dashed red arrows. The losses are high because the streamlines of the transverse component are refracted at the hollow core boundary and further in the empty space are directed radially. (right) Negative curvature hollow core micro - structured fiber. The streamlines of the transverse component are deviated by the optical vortices in the cladding, where the transverse component of the Poynting vector turns into zero, which leads to loss reduction

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.

Fig. 17

The map of the radial component of the Poynting vector for all solid band gap fibers with one row of rods (Λ=12μm). (left) λ=1μm ; (right) λ=1.5μm

Fig. 18

The streamlines of the transverse component of the Poynting vector (grey lines) and the lines of zero intensity ${\vec {P_{t}}= (P_{x} = 0 (red lines) ; P_{y} = 0 (green lines))}$ for all solid band gap fiber with a cross section shown in Fig. 17 at a wavelength of λ=1μm.

Fig. 19

The streamlines of the transverse component of the Poynting vector (grey lines) and the lines of zero intensity ${\vec {P_{t}}= (P_{x} = 0 (red lines) ; P_{y} = 0 (green lines))}$ for all solid band gap fiber with a cross section shown in Fig. 17 at a wavelength of λ=1.5μm.

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 micro - structured 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 power is possible in a certain region of the parameter space.

• Coherent structures in the continuous-discrete systems, such as, e.g., arrays of waveguides with disperive noninear media, can form light bullets - waves localised both in space and time. These structures can be rotating - vortex light bullets.

• Different structures of vortices arising in the Poynting vector transverse component of the core modes of micro - structured optical fibers can lead to 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.

Appendix

Numerical methods for multicore fibers

Compact implicit finite-difference scheme

We consider the periodic boundary value problem for a system of linearly coupled NLSE

$$ i\frac{\partial U_{k}}{\partial z} = \frac{\beta_{2}}{2}\frac{\partial^{2} U_{k}}{\partial t^{2}} - f_{k}, \ f_{k} = \gamma |U_{k}|^{2} U_{k} + \sum_{m=1}^{N} C_{k,m} U_{m}, \\ $$

(38)

$$ U_{k}(0,t) = U_{0,k}(t), $$

(39)

$$ \ U_{k}(z,T) = U_{k}(z,-T), $$

(40)

for complex functions U_k=U_k(z,t), k=1,...,N, where z∈[0,L] is the spatial variable and t∈[−T,T] is temporal variable. The functions U_k represent the envelopes of the pulses propagating along a multi-core fiber with N cores.

In the paper [159], the compact finite-difference scheme of order O(h²+τ⁴) proposed in [47, 48] was modified for a system of linearly coupled NLSE (38)

$$\begin{array}{*{20}l} i\frac{U^{n+1}_{k} - U^{n}_{k}}{h}=\frac{\beta_{2,k}}{2}\Lambda\left(\alpha U^{n+1}_{k} + (1-\alpha)U^{n}_{k}\right) \\- \left(\frac{f_{k}^{n+1} + f_{k}^{n}}2 + \frac{\tau^{2}}{12}\Lambda f_{k}^{n} \right), \ k = 1,...,N. \end{array} $$

(41)

In this case the right-hand side f_k of each equation contains the nonlinear cubic term and linear couplings and is written as

$$ f^{n}_{k} = \gamma |U^{n}_{k}|^{2} U^{n}_{k} + \sum_{m=1}^{N} C_{k,m}U^{n}_{m}. $$

(42)

Here β_2,k=β₂ is the group-velocity dispersion parameter for the mode k, γ_k=γ is the Kerr parameter, C_k,m=C is the coupling parameter.

Here the weight

$$\alpha=1/2-i/(3r)+c h, r=2\beta_{2}h/\tau^{2}, c>0$$

to obtain strong stability and the second order differential operator in time is defined in the standard way, so

$$\Lambda U^{n}_{k,j}=\frac{U^{n}_{k,j+1}-2U^{n}_{k,j}+U^{n}_{k,j-1}}{\tau^{2}},$$

where z_n=nh,t_j=jτ. It can be shown that the small change in the weight α is equivalent to the addition of an artificial second-order dissipative term ch²∂³U_k/∂z∂t² to Eq. (41). The analytic proof of an approximation order O(h²+τ⁴) of the scheme (41) is presented in [159].

Equation (41) is an implicit nonlinear scheme, and must be solved iteratively. To obtain the iterations, a simple method of linearization, similar to the two-step procedure used in the predictor-corrector schemes, is done for each fixed layer n on an evolutional variable z. Specifically, we define the sequence V₀,V₁,⋯,V_s (s is the iteration number), that converges to the solution Uⁿ⁺¹ 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]

$$\begin{array}{*{20}l} &\left(iI-\frac {\beta_{2} \alpha h}{2} \Lambda \right)\frac{V_{k,s+1} - U^{n}_{k}}{h} - \frac{\beta_{2}}{2}\Lambda U^{n}_{k} + \\ &+\left[\frac{\gamma |V_{k,s+1}|^{2} V_{k,s+1} + \sum_{m=1}^{N} C_{k,m} V_{k,s+1} + f_{k}^{n}}2 \right. \\ &+ \left. \frac{\tau^{2}}{12}\Lambda f_{k}^{n} \right] = 0 \end{array} $$

(43)

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⁻⁸ 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:

$$ \frac{\partial U}{\partial z} = (\hat{D} + \hat{N}) U, $$

(44)

where U=(U₁,...,U_N)^T is a vector of complex variables, $\hat {D}$ is a linear differential operator defining dispersion and linear couplings among variables U_k, k=1...N, and $\hat {N}$ is a nonlinear operator. These operators have the following structure for the system of equations (38) considered above:

$$ \hat{D} = i\left[C - I\frac{\beta_{2}}{2} \frac{\partial^{2}}{\partial t^{2}}\right], $$

(45)

$$ \hat{N} = {diag}\left\{i\gamma|U_{1}|^{2}, i\gamma|U_{2}|^{2},...i\gamma|U_{N}|^{2}\right\}, $$

(46)

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 operator $\hat {D} = 0$. In the second step, only linear effects are taken into account, and therefore the nonlinear operator $\hat {N} = 0$. In practice, the following symmetrized form of the algorithm is used to achieve the second order approximation:

$$ {}U(z + h, t) \approx \exp\left(\frac{1}{2}h\hat{N}\right) \exp(h\hat{D}) \exp\left(\frac{1}{2}h\hat{N}\right) U(z, t), $$

(47)

which can be represented as follows in order to obtain the final solution at a distance L=hM_z:

$$\begin{array}{@{}rcl@{}} U(L, t) \approx \exp\left(-\frac{1}{2}h\hat{N}\right) \left(\prod^{M_{z}}_{m=1} \exp(h\hat{N}) \exp(h\hat{D}) \right) \\ \exp\left(\frac{1}{2}h\hat{N}\right) U(0, t). \end{array} $$

(48)

When SSFM is generalized on the NLSE system (44), the calculation of the nonlinear diagonal operator $\exp (\alpha h\hat {N})$, α={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 operator $\hat {D}$ is a non diagonal matrix, therefore, calculation of the matrix exponential $\exp (h\hat {D})$ is needed for implementation of linear step. The following scheme was used to calculate it. The matrix exponential $\exp (h\hat {D})$ is calculated in Fourier space, using the representation

$$ \exp(h\hat{D}) \tilde{U}(z, t) = F^{-1}_{t} \exp[h\hat{D}(-i\omega)] F_{t} \tilde{U}(z, t), $$

(49)

where F_t denotes the Fourier transform operation, $\tilde {U}(z, t)$ is a vector with intermediate values of N variables that can be written as

$$ {}\tilde{U}(z, t) = \left(\prod^{n}_{m=1} \exp(h\hat{N}) \exp(h\hat{D}) \right) \exp\left(\frac{1}{2}h\hat{N}\right) U(0, t), $$

(50)

where z=nh, and the operator $\hat {D}(-i\omega)$ was obtained from the dispersion operator (45) by replacing derivative ∂/∂t by −iω, where ω is the frequency in Fourier space. The resulting matrix $\exp [h\hat {D}(-i\omega)]$ 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

$$\begin{array}{@{}rcl@{}} r_{l,m}(X) = \left[\sum\limits_{j=0}^{m} \frac{(l + m - j)! m!}{(l + m)!(m - j)!} \frac{(-X)^{j}}{j!} \right]^{-1} \\ \left[\sum\limits_{j=0}^{l} \frac{(l + m - j)! l!}{(l + m)!(l-j)!} \frac{X^{j}}{j!} \right] \approx \exp[X], \end{array} $$

(51)

where X is an N×N matrix, and represents the rational polynomial with the numerator p_l,m(X) of degree l and the denominator q_l,m(X) of degree m, providing the order of accuracy O(X^l+m+1)=e^X−r_l,m(X). The diagonal approximants r_m,m(X)=r_m(X)=p_m(X)/q_m(X) are preferred in practice, because r_l,m with l≠m is less accurate than r_j(X) with j= max(l,m), though r_j(X) can be evaluated with the same cost. In addition, the spectral radius ρ of approximant r_m(X) is less than 1, if the eigenvalues of the matrix X lie in the open left complex half-plane [165]. This property is important for the numerical solution of differential equations.

It is convenient to use the relation e^X=(e^X/σ)^σ in the calculation of the matrix exponential (scaling and squaring method), where σ is an integral power of 2, σ=2^s. This approach can significantly increase the accuracy of computing the matrix exponential e^X in comparison with the direct application of the numerical approximation formula (51) for matrices with large norm [164]. First, we must choose the parameter σ such that the norm of the matrix X/σ has become the order of 1. Then the Padé approximation is calculated (51) for the normalized matrix, which is raised to the power σ.

Higham proposed the scaling and squaring method based on a novel error analysis [122]. This method is highly accurate and fast due to an adaptive choice of Padé approximant degree depending on the matrix norm (approximants of degrees m=3, 5, 7, 9, 13 are used), a smaller number of scaling and squaring operations and matrix multiplications. We implemented this algorithm for computing the matrix exponential in our modification of the SSFM.