Vortical light bullets in second-harmonic-generating media supported by a trapping potential

We introduce a three-dimensional (3D) model of optical media with the quadratic ($\chi ^{(2)}$) nonlinearity and an effective 2D isotropic harmonic-oscillator (HO) potential. While it is well known that 3D \chi^2 solitons with embedded vorticity ("vortical light bullets") are unstable in the free space, we demonstrate that they have a broad stability region in the present model, being supported by the HO potential against the splitting instability. The shape of the vortical solitons may be accurately predicted by the variational approximation (VA). They exist above a threshold value of the total energy (norm) and below another critical value, which determines a stability boundary. The existence threshold vanishes is a part of the parameter space, depending on the mismatch parameter, which is explained by means of the comparison with the 2D counterpart of the system. Above the stability boundary, the vortex features shape oscillations, periodically breaking its axisymmetric form and restoring it. Collisions between vortices moving in the longitudinal direction are studied too. The collision is strongly inelastic at relatively small values of the velocities, breaking the two colliding vortices into three, with the same vorticity. The results suggest a possibility of the creation of stable 3D optical solitons with the intrinsic vorticity.


I. INTRODUCTION
The second-harmonic-generating (χ (2) , i.e., quadratic) nonlinearity is one of basic features of photonic media. Among other fundamental effects, the χ (2) interactions give rise to two-color solitons of diverse types [1][2][3][4], [5]. In particular, the use of this nonlinearity is critically important for the creation of stable two-and three-dimensional (2D and 3D) solitons, because, in contrast to the cubic (χ (3) ) self-focusing, the χ (2) interaction between the fundamental-frequency (FF) and second-harmonic (SH) fields does not induce the wave collapse [6], which destabilizes multidimensional solitons in χ (3) settings [4,7]. Due to this circumstance, the χ (2) solitons were originally created, in the spatial domain, as stable (2+1)-dimensional beams [8]. The possibility of the creation of fully localized 3D spatiotemporal solitons is suggested too by the absence of the collapse [9][10][11][12]. This has not been achieved yet in experiments, the best result being a spatiotemporal soliton which is self-trapped in the longitudinal and one transverse directions, while the confinement along the other transverse coordinate was provided by the waveguiding structure [13,14]. The χ (2) nonlinearity is also promising for the creation of ultrashort few-cycle optical solitons [15], The creation of self-trapped beams with intrinsic vorticity is another natural possibility in the 2D setting, with the SH and FF fields carrying topological charges m and m/2, respectively. For even m, these modes are classified as solitary vortices with topological charge m/2. States with odd values of m are not possible, as the topological charge of the FF component cannot take half-integer values (vortices with a half-integer optical angular momentum, in the form of mixed screw-edge dislocations, can be created by passing the holding beam through a spiral-phase plate displaced off the beam's axis [16]). Unlike the fundamental χ (2) solitons with m = 0, solitary vortices in the free space are always unstable against azimuthal perturbations, which split them into sets of separating segments, as shown both theoretically [17][18][19][20] and experimentally [21]. A similar splitting instability of vortices was predicted in the framework of the three-wave (alias Type-II) χ (2) system, which includes two components of the FF field [22,23].
Vortical solitons are well known too as solutions to the 2D and 3D nonlinear Schrödinger (NLS) equation with the cubic (χ (3) ) self-focusing nonlinearity [24]. They are also unstable against azimuthal perturbations, which develops faster than the above-mentioned collapse, i.e., the vortex splits into fragments, which later suffer the intrinsic collapse [4]. A method for the stabilization of the vortices was elaborated in the framework of the 2D and 3D Gross-Pitaevskii (GP) equations for atomic Bose-Einstein condensates (BECs) with attractive inter-atomic interactions, which are similar to the NLS equations for photonic media with the self-focusing χ (3) nonlinearity [25]. It has been shown in detail theoretically [26][27][28][29][30][31][32] that the axisymmetric harmonic-oscillator (HO) potential stabilizes the fundamental solitons in their whole existence domain, and vortices with topological charge 1 in a part of the region where they exist.
In optical media, the effective 2D trapping potential can be realized too, in the form of an appropriate transverse profile of the local refractive index [5]. The stabilization of 2D vortex solitons by means of the isotropic HO potential in the medium with the χ (2) nonlinearity was reported in our recent work [33]. A 2D photonic crystal can be used to build this setting, with the χ (2) nonlinearity provided by the poled liquid [34] or solid [35] material filling its voids. In fact, the profile of the effective potential does not need to be exactly parabolic (HO), due to the strong localization of the trapped modes. Our recent analysis [33] demonstrates the stabilization of both two-color χ (2) vortices and single-color ones, only the SH component being present in the latter case. Stability regions were found for both types of the modes. In particular, while the single-color states are actually linear modes supported by the trapping potential, their stability is determined by the χ (2) interactions with small perturbations in the FF field. Two-color trapped vortices exist above a finite threshold values, which is identical to the instability boundary of the single-color vortex mode.
A natural extension of the previous analysis, which is the subject of the present work, is to consider a possibility of the stabilization of two-color 3D vortex solitons ("vortical light bullets" [3]) by the 2D axisymmetric trapping potential, which may open the way to creating such self-trapped spatiotemporal modes. Thus far, they evaded attempts of the experimental observation (it is obvious that, unlike the 2D setting, single-color localized modes cannot exist in the 3D χ (2) system with the 2D potential). The 3D model, based on the system of coupled equations for the FF and SH fields, is introduced in Section II. In particular, an essential parameter, specific to the 3D model, is the ratio of the group-velocity-dispersion (GVD) coefficients for the FF and SH waves. In fact, essentially the same system of GP equations for atomic and molecular wave functions, applies to the BEC in atomic-molecular mixtures, where the role of the GVD ratio is played by the ratio of the atomic and molecular masses [36][37][38][39], hence the mechanism of the stabilization of the 3D vortex solitons trapped in the 2D HO potential may be implemented in the BEC mixture too. In Section II we also formulate the variational approximation (VA) for stationary states.
Stationary solutions for the trapped 3D vortical modes are considered in Section III. The results are obtained by means of numerical methods, and in a semi-analytical form based on the VA. In those cases when the underlying ansatz is appropriate, the accuracy of the VA is very good, in comparison with numerical results. Depending on the mismatch parameter, the vortex modes may exist above a finite threshold value of the norm (total energy), which can be explained in a simple form. The stability of the vortex solitons is considered, by means of systematic simulations of the perturbed evolution, in Section IV. The vortex modes are stable below a certain critical value of the norm. Above the instability boundary, the vortex develop oscillations, maintaining its vorticity and demonstrating periodic breaking and recovery of the axial symmetry of its shape. Collisions between vortex solitons moving in the longitudinal direction are also considered in Section IV, a noteworthy result being that, at sufficiently small velocities, two colliding vortices may break up into three. The paper is concluded by Section V.

II. THE MODEL AND THE VARIATIONAL APPROXIMATION
The equations for the evolution of the FF and SH field amplitudes, u (x, y, z, τ ) and v (x, y, z, τ ), in the 3D setting, with propagation distance z, transverse coordinates (x, y), and reduced time τ ≡ t − z/V are taken in the known form [1][2][3] (V is the common group velocity of the FF and SH waves, assuming that the group-velocity mismatch between them may be made equal to zero by means of an appropriate Galilean boost): where the asterisk stands for the complex conjugate, q is the real mismatch coefficient, the anomalous-GVD coefficient at the FF is normalized to be 1, and D is the ratio of the SH and FF GVD coefficients. In the case of D = 1, the diffraction-dispersion operators for the FF and SH components feature the spatiotemporal isotropy [9]. As said above, the axisymmetric modulation of the refractive index is represented by a spatially-isotropic HO potential, U (x, y) = Ω 2 /2 x 2 + y 2 . By means of a further rescaling, we fix Ω = 1/2, while q is kept as a free constant. Stationary modes are characterized by their total energy (norm), which is dynamical invariant of the system. Three other invariants are the angular momentum: the linear momentum in the τ direction: and the Hamiltonian: Equations (1) can be derived from the corresponding action, S = Ldz, with Lagrangian This fact may be used to develop the VA for spatiotemporal vortical solitons with propagation constant k, based on the following natural ansatz, which corresponds to topological charge 1 and is written in terms of polar coordinates (r, θ) in the plane of (x, y): Here (u 0 , v 0 ) and α −1/2 1,2 are, respectively, the amplitudes and radial widths of the two components, and β −1 is their common temporal width. The norm (2) of the ansatz is The substitution of ansatz (7) into the Lagrangian yields Using the relation following from Eq. (9) is cast into the form in which N 0 may be considered as a given constant, while u 0 , α 1 , α 2 , and β are treated as variational parameters: The system of variational equations produced by the latter Lagrangian, ∂L/∂u 0 = ∂L/∂α 1 = ∂L/∂α 2 = ∂L/∂β = 0, take a cumbersome form, which, nevertheless, admits a numerical solution: A characteristic example is presented in Fig. 1, which shows a plot of |u(x, y, τ )| in cross sections drawn through τ = L τ /2 (a) and y = L/2 (b) at N = 40, q = 1, D = 1, Ω = 0.5. In particular, Fig. 1(a) shows the vortex structure per se, with amplitude |u| vanishing at the center, and Fig. 1(b) shows the self-trapped structure in the τ direction. For the same parameters, the numerical solution of the variational equations (11) yields u 0 = 0.704, α 1 = 0.268, α 2 = 0.507, and β = 0.438. To compare the full numerical results with the prediction of the VA, the solid curve in Fig. 2(a) shows the numerically found profile of |u(x − L/2, y, τ )| at y = L/2, τ = L τ /2, and the dashed curve shows the VA-predicted counterpart of this profile, taken as per Eq. (7), i.e., u 0 |x| exp −α 1 x 2 . Further, the solid curve in Fig. 2(b) shows the temporal profile, |u(x, y, τ − L τ /2)|, at fixed x = 1.663 and y = L/2, and the dashed curve shows its VA-predicted counterpart, u 0 |x| exp −α 1 x 2 / cosh(βτ ) at the same point, x = 1.663. It is seen that the VA predicts the numerical results in a practically exact form.
Systematic results obtained for the vortex-soliton family in the direct numerical and VA forms are summarized in Fig. 3, in which panel (a) shows the slope of the vortical profile at its center, slope ≡ r −1 |u(r, τ )| at r → 0, as a function of the total energy (norm) N , for two values of the mismatch, q = 1 and 0.6. In the latter case, the vortex soliton disappears discontinuously at a critical value N = N c1 ≈ 19.5, as demonstrated by the numerical results and the VA alike, and the vortices do not exist at N < N c1 . Further, Fig. 3(b) shows the critical value as a function of the mismatch. Note that N c1 vanishes at q ≥ 1. It is relevant to note that a finite threshold value of the norm, N = N (2D) c1 , necessary for the existence of two-color vortical solitons in the 2D version of Eqs. (1), was found in our recent work [33]. As well as in the present model, N (2D) c1 vanishes at q ≥ 1. The difference from the 3D case is that stationary 2D single-color vortices trapped in the HO potential, with the zero FF component (i.e., these are linear SH modes), exist at N < N c1 . Obviously, this is not possible in the 3D case, as the 2D potential cannot trap 3D linear modes. The vanishing of N (2D) c1 at q ≥ 1 was explained analytically in our recent work [33]. In fact, the same explanation applies here, as the problem of the presence or absence of the finite threshold reduces to the existence or nonexistence of the solitons with a vanishingly small norm, N → 0. Obviously, in this limit the radial width of the mode trapped by the HO potential remains finite in the (x, y) plane, while the temporal width (in the τ direction) becomes indefinitely large. Thus, the 3D problem for N → 0 actually reduces to its 2D counterpart, for which an explanation of the absence of the threshold at q ≥ 1 was recently given in work [33].
We have also studied vortex solitons in the spatiotemporally anisotropic system, with D < 1 in Eqs. (1) (usually, the anomalous GVD is weaker at the SH [4], therefore it makes sense to consider the case of D < 1). In particular, Fig.  3(c) shows, for this case, the vortex' slope at the center, defined in Eq. (13), as a function of N at q = 1 and 0.6 for the limit case of D = 0 (zero GVD at the SH component). Note that the vortex soliton disappears continuously both for q = 0.6 and 1 at D = 0 [i.e., N c1 (q = 0.6, D = 0) = 0], but it disappears by a jump at q = 0, i.e., N c1 (q = 0, D = 0) is finite.
The VA approximation poorly agrees with the numerical results for D = 0, which is explained by the fact that ansatz (7), postulating equal temporal widths of the FF and SH components of the vortex soliton, is not appropriate in this case. For D < 0, vortex-soliton solutions could not be found, in agreement with the known general principle that they do not exist in this case [4,41].

IV. STABILITY AND DYNAMICS OF THE VORTEX SOLITONS
Systematic simulations demonstrate that the vortex solitons found above are stable in a large part of their existence domain. However, they are subject to an oscillatory instability at N > N c2 , when the total energy (norm) of the soliton exceeds the respective critical value, N c2 . An example of such an instability is presented in Fig. 4(a), which shows the evolution of two azimuthal Fourier amplitudes, which are defined as at q = 1, D = 1 and N = 185 (obviously, for the stationary vortex u +1 = const, and u −1 = 0). In the course of the unstable evolution, amplitude u −1 (z) increases from zero, exhibiting oscillatory dynamics. Further, Fig. 4(b) shows the perturbed evolution of the same amplitudes for q = 1, N = 185, and D = 0 (zero GVD at the SH component). The difference from the picture shown in 4(a) is that, at D = 0, the oscillations are apparently chaotic. Figure 5 demonstrates that, in terms of the shape of the same vortex soliton whose evolution is illustrated by Fig.  4(a), the oscillatory perturbations induce an azimuthal modulational instability of the axisymmetric shape. The axial symmetry is partly destroyed and then restored periodically, as one can see from the comparison of Figs. 5 (a) and (b). In the former panel, the modulation depth attains its maximum exactly at point z = 110, where u −1 (z) has a maximum in Fig. 4 (a), and the axial symmetry is restored at z = 150, where u −1 (z) vanishes in Fig. 4(a). The period of the periodic breakup and retrieval of the axial symmetry is Z ≈ 80, the same as observed in the oscillations of u −1 (z) in Fig. 4(a). This dynamical regime resembles the ones featuring periodic splittings and recombinations of 2D vortices, which are trapped in the isotropic HO potential in the systems with χ (2) [33] and self-attractive χ (3) [29] nonlinearities. The critical value for the onset of the oscillatory instability, found from the direct simulations, is shown in Fig. 6(a) for several values of the mismatch, q = 0.2, 0.6, 1.0, 1.4, and 1.8, in the systems with the spatiotemporally isotropic diffraction-dispersion (D = 1), and with zero GVD of the SH component (D = 0). Stable vortex solitons exists in the broad region of N c1 < N < N c2 [recall N c1 is the existence threshold for the solitons, see Fig. 3(b)]. It is worthy to note that the stability region is somewhat larger for D = 0 than for D = 1.
The obvious Galilean invariance of Eqs. (1) with D = 1 (the case of the spatiotemporally isotropic diffractiondispersion operators) in the longitudinal (τ ) direction suggests a possibility to study collisions between moving solitons, if they are set in motion by means of the application of properly matched kicks to the FF and SH components, (u, v) → ue −iωτ , ve −2iωτ , where ω is an arbitrary boost parameter. Figure 7(a) shows a typical example of the simulated head-on collision of two vortex solitons with equal vorticities 1, kicked by ω = ±0.5. The collision is strongly inelastic, transforming the two colliding vortices into three, one quiescent and two symmetrically moving ones. We stress that each soliton produced by this inelastic interaction is a vortex with the same vorticity, 1, which was embedded into the initial modes, the total angular momentum (3) and the linear momenetum (4) being conserved. Inelastic collisions between vortices do not generate additional zero-vorticity solitons. The increase of the kick to ω = ±1 results in a quasi-elastic passage of the colliding solitons through each other, see Fig. 7(b). Figures  7(c) and (d) demonstrate the collision of two solitons with opposite vorticities, +1 and −1 kicked by (c) ω = ±0.5 and (d) ω = ±1.@The collision is inelastic for ω = ±0.5 and quasi-elastic for ω = ±1.

V. CONCLUSION
The objective of this work is to propose the model of optical media with the χ (2) nonlinearity, in which 3D vortical solitons, which are always unstable against splitting in the free space, may be stabilized by means of the effective 2D isotropic HO (harmonic-oscillator) potential, induced by an appropriate spatial modulation of the refractive index. A broad stability region for the 3D vortices has been found, with their shape accurately predicted by the VA (variational approximation). The vortices exist with the norm (total energy) exceeding a finite threshold, which vanishes at values of the mismatch q ≥ 1. The latter fact was explained by means of the comparison with similar results recently reported in the 2D version of the system in our work [33]. The stability region is limited by another critical value of the norm, which is much higher than the existence threshold. Above the stability border, the vortex oscillates, periodically breaking and restoring its axisymmetric shape. Collisions between vortices, which may move freely in the longitudinal (temporal) direction, were also studied. If the collision velocity is small enough, the collision transforms the two colliding vortices into three, with the same vorticity. The results reported in this work suggest that 3D optical vortex solitons may be created in the stable form.
The analysis can be naturally extended in other directions. In particular, it may be interesting to study the system with the Type-II χ (2) nonlinearity, i.e., the three-wave system with two FF components determined by their mutually orthogonal polarizations [2,3]. Further, it is relevant to consider a possibility to stabilize the 3D vortex solitons using a periodic (lattice) potential, instead of the HO one, cf. the analysis performed previously for the χ (2) [42] and χ (3) [43][44][45] nonlinearities. Also relevant may be the analysis of the stability of trapped vortices in a system including competing χ (2) and χ (3) nonlinearities, cf. the consideration of the free-space model in work [46].