Finding spatiotemporal light bullets in multicore and multimode fibers.

A two-level iterative algorithm for finding stationary solutions of coupled nonlinear Schrödinger equations describing the propagation dynamics of an electromagnetic pulse in multimode and multicore optical fibers of various structures was developed and tested. Using as an example the proposed analytical soliton solution which is localized in space and time, test calculations were performed, and the convergence of the algorithm was demonstrated.


Introduction
Multicore fibers (MCFs) and multimode fibers (MMFs) are currently used in a large number of scientific and technical applications. Despite the fact that initially they emerged as solutions to increase the throughput of fiber-optic communication lines, nowadays MCFs and MMFs are used in a wide range of photonics applications. MMFs have one core supporting a few number of optical modes, while MCFs consist of several small cores located under a common cladding (see Fig. 1). Such structures can be considered as nonlinear discrete physical systems, that are important both for basic research and for practical use. For example, it has recently been shown that propagation of high-power optical pulses through a MCF can lead to energy concentration in certain cores, leading to the coherent combination, amplification and compression of light pulses [1]. This nonlinear effect arising in MCFs can be used to create photonic devices, such as optical switches, as well as to develop new types of high-power fiber lasers. In addition, multimode and multicore fibers have great potential for applying to the effective transmission of information over long distances. Their use can serve as an alternative to the parallel laying of fiber cables, since the bandwidth in such fibers increases in proportion to the number of modes or cores [2,3]. Record performance for data transfer have recently been achieved with the help of a combined MCF-MMF approach. For example, in [4] using a 19-core fiber and 6 modes in each core, a record bandwidth of 10.16 petabits per second was achieved. Data transmission was carried out at a distance of 11.3 km.
Different spatio-temporal effects can be observed in MMFs. Recently, the effect of spatial beam self-cleaning was discovered [5,6]: above a certain threshold power, a cleaned powerful pulse can be obtained from an initial random optical field introduced in the MMF core. The Gross-Pitaevskii equation (GPE) is a suitable mathematical model for numerically describing such type of dynamics, which however is too complex for qualitatively describing this effect. It is possible to construct a simpler, low-dimensional model based on mode decomposition, which can provide some insight into the self-cleaning phenomenon (see [5]). As we shall see, a simple reduced "mean field" model can be suitable for the qualitative description of the steady-state generation in MMFs.
An interesting phenomenon associated with multicore fibers is the generation of light bulletshigh-energy pulses localized both in time and space, and propagating without distorting their shape along the fiber. Light bullets are formed due to the mutual compensation of diffraction, anomalous dispersion, and nonlinearity; they were recently observed in experiments [7][8][9]. In bulk media, small fluctuations in the intensity or width of the light beam can lead to collapse; in such case, light bullets are unstable. However, in discrete media, such as for example in MCFs, a similar effect is not observed [10]. In recent years, due to the increased interest in multicore light guides, many papers have appeared which are devoted to theoretical studies of light bullets. The propagation of the electromagnetic field in multicore and multimode fibers can be described by a system of coupled nonlinear Schrödinger equations (CNLSEs) or generalized Manakov equations [1,11]. Approximate analytical solutions in the form of space-time discrete solitons [12] were obtained for models of centrally symmetric multicore optical fibers based on the CNLSEs.
But in most cases, the search for stationary solutions of pulse propagation in MCFs and MMFs is only possible with the help of numerical methods. Among the methods for finding stationary solutions, the most common are: the Petviashvili method [13], as well as the imaginary time evolution method (ITEM) [14]. In the ITEM method, a solution of predetermined power is sought. Time t is replaced by the imaginary time it, for constructing an iterative process in the norm of L 2 . The disadvantage of this method is the inability to fix the propagation constant, which needs to be calculated at each iteration step. The Petviashvili stabilizing correction method, on the other hand, is suitable for finding solutions with a predetermined propagation constant, but the power of the obtained solution is not fixed. In the case of nonlinear photonic structures (for example, MCFs), the propagation constant parameterizes well the localized solution [13,15]. However, this algorithm is only applicable to MCFs with a regular structure and a large number of cores or for ring MCFs. In such cases one can reduce the relevant system to an uncoupled form by applying discrete Fourier transform by the discrete core space. For more complex MCFs, both of the above methods converge either to a trivial solution or to a distributed, non-localized solution.
Here we consider the case of a MCF with a centrally symmetric arrangement of the cores (see Fig. 1(c)). Since the central core introduces an irregularity into the structure, the use of standard iterative methods for finding a stationary solution in the form of a light bullet of nonzero energy is difficult in this situation.
In this work, we introduce a new two-level iterative algorithm, which allows for finding a pulsed stationary solution of the CNLSEs with different types of couplings among the cores of the MCF. The propagation of the found solutions was simulated, in order to demonstrate truly stationary dynamics. Next, we derive a new approximate analytical solution for ring MCFs with central core, and demonstrate the stability of such a solution. In addition, the analytical solution can also be used as initial distribution in order to start the iterative algorithm. Finally, we briefly investigate spatiotemporal bullet solutions for MMFs.

Mathematical model for MCF
The propagation of the electromagnetic field in a MCF can be described by using a system of CNLSEs, defined in the slowly varying envelope approximation. We write a system of (N C + 1) equations, which characterize the field propagation in a MCF with a centrally symmetric structure, as follows where A n is the envelope of the field propagating along the nth core, D n > 0 is the dispersion parameter of the nth core (we consider only the case of anomalous dispersion), γ n is the Kerr nonlinear coefficient, z is an evolutionary variable, t is a time, and N c is the number of cores on the periphery of the MCF. Further, for each core, we will only take into account the coupling with nearest cores: for example, for the cores at the periphery, the right and left neighbors, as well as the central one, will contribute to the evolution of the field. The coupling decreases in inverse proportion to the square of the distance between the cores, so this simplification is acceptable. The core located in the center is equally spaced from cores on the periphery, so that all members of the sum C n,m A m are taken into account. Due to the symmetry of the waveguide structure, the stationary pulse solution will also have a central symmetry: in the cores at the periphery, pulses will be characterized by one set of parameters (A 1 , γ 1 , D 1 ). The solution in the central core will be described by the set (A 0 , γ 0 , D 0 ). Therefore, we will only consider 2 equations, one for the central core, and one for a generic core on the periphery: Let us denote C 0,0 = β 0 , C 0,n = C 0 , C n,n+1 = C 1 , C n,n = β n , and introduce the dimensionless variables [16,17]: The final system of equations is written in the form (omitting primes):

Iterative method
In order to construct an iterative method, we will perform several transformations of the system Eq. (3). First, we will search for a stationary solution in the form u n (z, t) = u n (t) exp (iλ 2 z). After substituting this type of solution into the system, we divide the first equation by 2N c γ 0 γ 1 , and the second by 2, so that all coefficients of nonlinear terms are equal to one. Next, we calculate the Fourier transform of both equations of the system, and obtain whereũ n =ũ n (ω) and N(ũ n ) = F(|u n | 2 u n ) is the Fourier transform of the nonlinear term. As one can see, system Eq. (4) has an asymmetry of linear couplings. Next, we rewrite the system Eq.
(4) in a matrix form Mũ =ṽ: The two-level iterative method for searching a stationary solution of system Eq. (4) is schematically presented in Fig. 2. The purpose of the outer loop is to determine the the non-zero total energy E of the stationary solution, while the inner loop is responsible for the construction of the solution shape. As initial non-zero conditions u (0) 0 (t) and u (0) 1 (t), we used the approximate solutions that we will derive in the next section. For the total energy of the initial approximation , we define the right E max = 10 2 E 0 and the left E min = 10 −2 E 0 energy interval boundaries, within which we will search the energy of a stationary solution by means of the bisection method. The condition for stopping the outer loop iterations is provided by the criterion where ε is a given accuracy, and k is the outer loop iteration number. In calculations we used the value ε = 10 −11 . The inner iteration loop is necessary for approximating the shape of the stationary solution with a target value of energy E M , and consist of the transformation Eq. (5) coupled with a renormalization of the resulting solution energy: Here the index l is an iteration number of the inner loop. The stopping criterion is given by the stabilisation of the total energy of the field M −1ṽ(k,l) (before renormalization). To calculate the nonlinear vector of the right-hand sideṽ (k,l) , it is necessary to perform the inverse Fourier transform after each inner iteration in order to calculate N(ũ (k,l−1) n ). In general, this slows down the algorithm, so, for example, when searching for solutions of uncoupled systems of equations, the Petviashvili method will have an advantage in terms of computational speed. However, the advantage of the present algorithm is its applicability to the presently considered type of irregular waveguide structure, as discussed in the introductory Sec. 1. The application of the method to an arbitrary system of equations should not be confusing, since the matrix M for any system with linear couplings is diagonal, and generalization of the algorithm for the systems with more than 2 equations is quite straightforward. The application of the algorithm, generalized to the case of many (i.e., at least a dozen) equations, made it possible to find stationary solutions in [12,18].

Analytic approximate stationary solution
For a verification of the numerical iterative algorithm, and in order to define applicable initial conditions, we shall derive now an analytical approximate solution. We search a solution in the form of a light bullet, where most of the energy is concentrated in the central core, and low energy pulses are distributed at the periphery. In addition, let us assume that a solution in the form of a light bullet corresponds to a spatiotemporal soliton, the envelope of which will be in the form of a hyperbolic secant. The solution for the central core is whereas for the peripheral core: where ρ k are the pulse widths in the central and peripheral cores, respectively, f k are test functions, and A, B are functions that approximate the amplitude of the stationary solution. The solution Eq. (8) is a fundamental soliton solution, and the solution Eq. (9) is a localized optical pulse [12]. The pulse width of the central core can be found from the condition of a mutual compensation of dispersive and nonlinear effects, specifically Since the pulse at the periphery is not a fundamental soliton, we will also find its width by the method of indeterminate coefficients. By substituting the solutions Eq. (8) and Eq. (9) into the investigated system of Eq. (3) at the point t = 0, and using the method of undetermined coefficients for polynomials of λ, we obtain the coefficients α k , β k and r k : 2 ln 2α 1 , (2) , where s = −κ 2 + κ 4 + 8α 2 1 . The disadvantage of proposed the analytic expressions is the need to fulfill the condition s > 1, which leads to the limitation on the core number N c < 8 γ 1 γ 0 − 1 . For the analytic approximate solutions of Eq. (3), we obtain the following expressions for the total energy and the Hamiltonian: where d 1 = 5.69398194001592, d 2 = 7.37436321044295.

Analysis of stationary solutions of the CNLSEs
We carried out a comparative analysis of our analytical approximate and numerical solutions of the CNLSEs Eq. (3). We considered the following parameter values: D 0 = D 1 = 1, γ 0 = 1, γ 1 = 2. Without any loss of generality, the value κ = 2/ √ N c was chosen. Moreover, we took the time interval t ∈ [−30; 30] uniformly divided by M = 2 14 nodes.
In Fig. 3, we illustrate the results of the comparison of the obtained pulse profiles in the central (a,c) and peripheral (b,c) cores, for different values of the parameter λ, and for the core numbers N c = 2 and N c = 6, respectively. In the central core, a good agreement between analytical and numerical solutions is observed, both in amplitude and in pulse width. The relative error of the discrepancy among approximate analytical and numerical solutions was ≈ 10 −3 . Note that the case with N c = 6 peripheral cores can be considered as a 7-core hexagonal MCF also. For the peripheral cores, a good agreement is only observed in amplitude, for small values of the parameter λ; whereas a slight discrepancy in the pulse widths is obtained. This can be explained by the use of different approaches for the calculation of the analytical solution. Specifically, the pulse width in the central core is uniquely calculated from the condition of mutual compensation of nonlinear and dispersive effects of the fiber Eq. (10). On the other hand, pulse width for the periphery cores is obtained by using the polynomial Eq. (9), the highest power of which is λ. Therefore, for small values of λ, possible discrepancies may arise between the widths of the analytical and numerical solutions. Whereas for large values of λ, the relative error decreases to an order similar to 10 −3 . The pulse amplitudes for both the central and peripheral cores are given by approximating linear combinations of inverse power functions with an accuracy of O λ −7 . Therefore one obtains a good agreement in both cases.
The relative error between analytical and numerical solutions for the pulses in the central and peripheral cores is presented in Fig. 4, as a function of the propagation constant λ. Also, the discrepancy between the total energy E and the Hamiltonian H is shown. Figure 5 presents an example of propagation of the analytic solution and of the numerical solution. The analytical solutions can give adequate approximation of numerical solutions in many cases. All approximation errors decrease with the growth of parameter λ, except for the relative error for the field in peripheral cores, which stabilizes to a certain level. This level increases with the growth of the total number of cores. For small values of the parameter λ, analytical expressions poorly describe the stationary solution; nevertheless they can still be used as an initial approximation for the proposed numerical algorithm. The discrepancy is caused by a change of soliton solution type for the MCFs: for small values of λ the stationary solution ceases to be an optical bullet, and most of the energy concentrates in the peripheral cores. This case will be considered later.  Let discuss now in more detail the properties of stationary solutions, which are numerically found by the proposed algorithm. As we have mentioned before, for each peripheral core number N c some value λ * = λ * (N c ) exists, for which the stationary solution with λ > λ * is a light bullet where most of its energy is concentrated in the central core. Whereas the stationary solution with λ < λ * represents a distributed solution. Figure 6 shows the total energy of the system E = E(u 0 ) + E(u 1 ) and the Hamiltonian H. Here these solution families are well separated by a discontinuity of the energy E = E(λ) (for N c = 2), or by a discontinuity of the energy derivative (for N c = 4 and higher values). To determine the stability of the found solution we can use the Vakhitov-Kolokolov criterion [19], which states that a positive slope linear dependence of the total energy E of the system on the parameter λ 2 indicates the stability of the solution. Also, it should be noted that in discrete systems a collapse does not occur for a negative value of the Hamiltonian H [20]. Thus, based on Fig. 6(a), we can state that all stationary solutions for N c = 2 are stable, and can propagate without a shape change. On the other hand, for N c ≥ 4, a value λ s = λ s (N c ) exist, such that the solution obtained for λ * < λ < λ s is unstable.   Fig. 7, we pointed out by stars on the total energy E = E(λ, N c ) curve some specific solutions and drew the corresponding pulse dynamics along fiber length z. Solutions 1) and 2) correspond to the distributed and soliton stationary regimes for a MCF with N c = 2 peripheral cores. Both solutions are stable due to the Vakhitov-Kolokolov criterion, and propagate without change of shape in numerical simulations performed by the split-step Fourier method adapted for CNLSEs [21]. Solutions 3) -5) belong to different ranges of the parameter λ. Regime 5) is a stable soliton solution, while regime 4) corresponds to an unstable soliton solution, since the Vakhitov-Kolokolov criterion is not met in this case. Nevertheless, the dynamics of regime 4) show that for a given value of λ there is a stable distributed solution that was not captured by the algorithm, due to insufficiently close initial data. Finally, solution 3) belongs to a family of stable distributed solutions.
In Fig. 8, we present the energy E 0 = E(u 0 ), the width (FWHM) w 0 , and the peak power P 0 for the central core only. Please note that the width of all solutions of the "light bullet" type is the same, for different values of the parameter λ and of the number of cores N c . In this case, the energy and the peak power of the solutions increase in proportion to each other. A different behavior is only observed for N = 2, where another type of solution exists.

Mathematical model for MMF
The proposed two-level algorithm can also be applied to systems describing optical pulse propagation in MMFs (see [5]): where ψ 1,2 = N 1,2 exp(−iθ 1,2 ), and N 1,2 = ψ 1,2 2 measures the number of photons in the fundamental mode and the generic higher-order-mode (HOM), respectively, K is a constant Fig. 9. Stationary solutions of Eq. (14) for different values of the parameter λ for the fundamental mode (left) and generic HOM mode (right), respectively. coupling coefficient originating from nonlinear four-wave mixing terms, and U 1,2 is a constant determined by overlap integrals. Also, similar systems appeared in problems connected with Bose-Einstein condensates [22] and nonlinear directional couplers [23]. Figure 9 presents the stationary solutions of Eq. (14) for λ = 4 and λ = 16, in the case with U 1 = 1, U 2 = 2, D 1 = D 2 = 1, K = 1. The intercomparison of the agreement between the solutions of Eq. (4) and numerical solutions of the GPE is left for a subsequent publication.

Conclusion
We derived and verified the validity of a numerical algorithm for finding a stationary solution of a system of nonlinear differential equations, describing the dynamics of an electromagnetic pulse in a multicore optical fiber with a centrally symmetric structure. The proposed two-level iterative method allows one to numerically find solutions for an irregular structure of the fiber cores, which is similar to light bullets observed in experiments [9]. The application of the algorithm presented in this paper made it possible to find stationary solutions in [12,18]. For the considered multicore fiber involving a radial structure with an irregularity in its center, a stationary soliton solution which is localized both in time and space was found. In addition, we confirmed that the found family of solutions is stable upon propagation.