Collisions of Light Bullets with Different Circular Polarizations

Collisions of left- and right-polarized spatiotemporal optical solitons have been numerically simulated for a locally isotropic focusing Kerr medium with anomalous chromatic dispersion. The stable propagation of such ``light bullets'' in a moderate nonlinear regime is ensured by a transverse parabolic profile of the refraction index in a multimode waveguide. The transverse motion of centers of mass of wave packets in such systems occurs on classical trajectories of a harmonic oscillator, whereas the motion in the longitudinal direction is uniform. Therefore, collisions of two solitons can be not only head-on but also tangential. An inelastic collision of two solitons with opposite circular polarizations can result either in two binary light bullets combining the left and right polarization or in more complex bound systems. DOI: 10.1134/S0021364024600691


Introduction
In nonlinear optics, spatiotemporal optical solitons have been studied for several decades (see [1][2][3][4][5][6][7][8][9] and numerous references therein).Such stably propagating coherent structures confined in three dimensions are called light bullets [1][2][3].Difficulty of this problem is that the nonlinear Schrödinger equation used as a basic mathematical model for a quasi-monochromatic light wave in a Kerr medium does not have stable soliton solutions in homogeneous three-dimensional space.A wave packet either spreads due to dispersion and diffraction or collapses due to nonlinearity (see [8][9][10][11] and references therein).For this reason, an additional stabilizing factor -saturated nonlinearity, spatial inhomogeneity, etc. -is necessary (see, e.g., [4,6,7,[12][13][14][15][16][17][18][19][20][21][22]).In particular, the effective transverse potential keeps a moderately nonlinear wave packet from diffractive (transverse) spread in a multimode waveguide with a smooth (approximately parabolic) profile of the refractive index, and anomalous chromatic dispersion in the longitudinal (temporal) direction is stably balanced by nonlinearity becoming effectively one-dimensional [2,3,5,8,9].In this case, the transverse degrees of freedom are not frozen, in contrast to the situation in single-mode optical fibers.The established balance results in the stability of light bullets with not too high total energy.Furthermore, the center of mass of such a soliton can be displaced in the transverse direction in a well-known way (on classical oscillator trajectories) so that its motion regime is intermediate between the rigid transverse fixation in a single-mode optical fiber and complete freedom in a homogeneous space.
However, the conventional scalar nonlinear Schrödinger equation describes only one polarization of light, linear or circular.A locally isotropic Kerr medium allows the nonlinear interaction between waves with two polarizations.In this case, the system of two incoher- * Electronic address: ruban@itp.ac.ru ently coupled nonlinear Schrödinger equations for slow complex amplitudes of left and right circularly polarized waves is appropriate [23].Correspondingly, it is necessary to consider combined solitons, i.e., stable nonlinear solitary structures with both circular polarizations generally in an arbitrary proportion.It is also important to know processes accompanying collisions of threedimensional optical solitons with different polarizations.To the best of my knowledge, such problems were studied only for one-dimensional systems.Very nontrivial results were obtained for inelastic collisions of one-dimensional vector solitons (see, e.g., [24][25][26] and references therein).Even more interesting dynamics should be expected in three dimensions.
A variational analysis performed in this work indicates the possibility of existence of moderately nonlinear binary light bullets.The three-dimensional numerical simulation confirms their stable nature.A number of numerical experiments on head-on and tangential collisions of three-dimensional light bullets with opposite circular polarizations are also carried out.Processes accompanying these inelastic collisions significantly differ from those occurring in the case of the interaction of identically polarized solitons.It is found that, depending on the parameters, the collision can result in two binary light bullets propagating at changed velocities away from each other and sometimes in the "merging" of colliding solitons, forming a more complex nonstationary structure.Wave collapse occurs at some initial parameters.

Main properties of the model
We consider a transparent optical medium with the focusing Kerr nonlinearity and the dispersion relation of linear waves k(ω) = ε(ω)ω/c with the anomalous dispersion k ′′ (ω) < 0 in a certain frequency range.In particular, these properties are inherent in many glasses.The known equation for the vector envelope of a weakly nonlinear quasimonochromatic light wave with the carrier frequency ω 0 in the paraxial approximation is taken as the main model.In the dimensional variables, this equation has the form where ζ is the coordinate along the beam, k ′ 0 = 1/v gr is the inverse group velocity of light in the medium, k ′′ 0 is the negative chromatic dispersion coefficient, ε(x, y, ζ) is the small inhomogeneity of the relative permittivity at the carrier frequency, and α(ω 0 ) and β(ω 0 ) are the positive coefficients of nonlinearity.Let τ = t − ζ/v gr be the "retarded" time.In terms of slow amplitudes ψ 1,2 (x, y, τ, ζ) of the left and right circular polarizations, the amplitude of the electric field is and light is then described by two coupled scalar nonlinear Schrödinger equations [23], similar to a binary Bose-Einstein condensate of cold atoms (with the change of variables ζ → t, τ → z).The appropriate rescaling leads to the dimensionless system where ∆ = ∂ 2 x + ∂ 2 y + ∂ 2 τ is the three-dimensional Laplace operator in the r = (x, y, τ ) "coordinate" space and g = 1 + 2β/α is the cross-phase modulation parameter, which is about 2 in the typical case of fast nonlinear response.It is worth emphasizing that the nonlinear interaction between two components is reduced to a simple noncoherent coupling through the coefficient g only in the basis in the form of circular polarizations.Such a coupling conserves the amount of each component N 1,2 = |ψ 1,2 | 2 dxdydτ and makes it possible to apply the highly efficient numerical split-step Fourier method (see below).Nonlinear terms corresponding to the socalled four-wave mixing would remain in the system with two linear or two elliptic polarizations.These terms are inconvenient both analytically (because they result in transitions between components) and numerically except for the obvious reduction ψ 2 = ψ 1 exp(iδ 0 ), which corresponds to the purely linear polarization.
The external parabolic potential U ∝ −ε(x, y, ζ) can be represented in the compact matrix-vector form The coupled nonlinear Schrödinger equations (3) constitute a hydrodynamic system for two interacting compressible perfect liquids whose flows are potential.The transition to the corresponding hydrodynamic variablesdensities I 1,2 and potentials of velocities ϕ 1,2 -is ensured by the Madelung transformation A similar system of equations but with defocusing nonlinearity was recently considered in [27].In the defocusing case, the positive "hydrodynamic" pressure supports a delocalized density background, and the main "soft" coherent structures are quantized vortices and domain walls between two circular polarizations.In this work, the case of focusing nonlinearity is considered, and of main interest are soliton-like objects where a negative hydrodynamic pressure (tension) is balanced by the "quantum" pressure (i.e., dispersion) in the longitudinal direction, whereas diffraction in the transverse direction slightly weakened by the hydrodynamic tension resists the transverse confining potential.
It is also important that the center of mass of any localized distribution of liquids in the quadratic external potential moves according to the equation of the classical oscillator (cf.[28]): where the double prime means the second derivative with respect to the variable ζ.In some cases, the trajectory of this motion can be quite nontrivial, e.g., in the rotating anisotropic potential when the transition to the rotating reference frame adds the Coriolis and centrifugal forces.The periodic dependence of the eigenvalues of the matrix Ŝ(ζ) providing conditions for the parametric resonance is also interesting.However, even a conventional anisotropic two-dimensional oscillator with the constant matrix Ŝ = Diag{κ 2 1 , κ 2 2 , 0} "yields" Lissajous figures.In this case, it is very substantial that the "internal" dynamics of the localized structure is "insensitive" to the motion of the center of mass.Indeed, let us consider a slightly "truncated" auxiliary system of equations without terms with F(ζ) and U 0 (ζ): It is easy to verify that the functions satisfy the complete system of equations (3) under the conditions The first two conditions are the equations of motion of the classical oscillator in the presence of the driving force F(ζ).Without loss of generality, it can be accepted that Ψ 1,2 correspond to the center of mass at rest; then, it becomes obvious that the motion of the center of mass does not affect the internal dynamics, as in the case of the conventional nonlinear Schrödinger equation [28].

Gaussian variational ansatz
This internal dynamics itself can be fairly complex, particularly, in the presence of both polarizations.In particular, a combined soliton including two components in the numbers N 1 and N 2 is certainly of interest.Such binary light bullets have not yet been studied.The properties of these objects in the simplest case of two symmetric wave packets with coinciding centers can be understood using the variational method.The system of equations ( 6) corresponds to the Lagrangian where the Hamiltonian functional has the form By analogy with the single-component nonlinear Schrödinger equation, trial functions in the form of Gaussian wave packets where Âs (ζ) and Bs (ζ) are unknown symmetric matrices and φ s (ζ) are unknown phases, are substituted into the Lagrangian.We emphasize that the matrices Âs (ζ) and Bs (ζ) are not generally assumed to be diagonal, in contrast to most of the works where the variational method is applied to the nonlinear Schrödinger equation.Gaussian integrals are easily calculated in the general form and lead to a Lagrangian system with a finite number of degrees of freedom: The interaction between the polarizations is obviously described by the last term.Naturally, the variables N s are integrals of motion.The equations for the matrices A s (ζ) and Bs (ζ) have the structure The explicit form of partial derivatives with respect to the elements of the matrix Â−1 s is easily determined using the relations d Tr M = −Tr( M 2 d M −1 ) and , where M is an arbitrary symmetric matrix.As a result, the system of ordinary differential equations is obtained in the compact matrix form In the general form, this system is still too complex for a detailed analysis.However, if only the diagonal matrices Ŝ(ζ), Âs (ζ), and Bs (ζ) are retained and Â−1 where slightly more convenient variables Ñs = N s / √ 8π 3 are introduced and the potential energy W has the form The corresponding dynamic system for the singlecomponent nonlinear Schrödinger equation is studied in detail (see, e.g., [2,3,8,9,29]; it is noteworthy that the longitudinal dependence in the diagonal variational ansatz is sometimes taken in the form 1/ cosh[τ /c(ζ)] rather than in the Gaussian form; the potentials W in this case are fairly similar).If κ 2 1 and κ 2 2 are independent of the evolution variable ζ, a local minimum of the function W determines the equilibrium configuration of a binary light bullet.In the axisymmetric case (where κ 2 1 = κ 2 2 = 1), it is obvious that a s = b s , and only four of six variables remain independent.In this case, the qualitative result is the same as for the single-component soliton: the function W has a stable equilibrium position in a certain region of not too large parameters Ñ1 and Ñ2 .If Ñ1 > Ñ2 , then a 1 > a 2 and c 1 > c 2 ; i.e., the weaker component is slightly more localized.This is natural because the nonlinear potential well (−I 2 − gI 1 ) for the weaker component at I 2 ≪ I 1 is deeper than the well (−I 1 − gI 2 ) for the stronger component due to the parameter g > 1.

Three-dimensional numerical simulation
To verify the conclusions provided by the variational analysis and to study collisions of solitons, the direct numerical simulation of the system of equations ( 3) by the standard second-order split-step Fourier method in the variable ζ was performed.The computational domain had sizes (4π) × (4π) in the transverse x and y directions and 8π or 12π in the longitudinal τ direction.The lattice step h = 4π/128 ∼ 0.1 in the (x, y, τ ) "space" together with the evolution step δζ = 0.002 ensured a sufficient high resolution of smooth wave fields (in the absence of collapse) and the conservation of the Hamiltonian with an accuracy of four to six decimal places in the entire propagation distance of several hundreds of dimensionless units of ζ.The experiment was stopped if a sharp and strong (several times) increase in the maximum intensity began shortly before the collapse although the spatial resolution in this case still remained acceptable.For this reason, the process of collapse was not simulated in this work.
Figures 1a-1c show the numerically obtained maximum intensities of the first and second components of the nonlinear binary wave packet with Gaussian initial conditions.As should be expected, detailed coincidence with the Gaussian approximation was not achieved.Degrees of freedom "disregarded" in the Gaussian were rapidly excited in computer experiments due to the imperfection of the initial conditions.However, it is important that these perturbations were not enhanced at long distances; consequently, the hypothesis of the stable propagation of binary light bullets was confirmed by the numerical test.
At periodic dependences of the coefficients of the quadratic potential such as κ 2 1 (ζ) = 1 + ǫ 1 cos(2ζ) and κ 2 2 (ζ) = 1 + ǫ 2 cos(2ζ) with small parameters ǫ 1 and ǫ 2 , the parametric resonance occurs in the system, as well as in the single-component case (see [30,31] and references therein).The computer simulation of Eqs.(3) indeed demonstrated the parametric "build-up" in the sizes of the binary soliton and a strong increase in its energy up to the situation where the transverse size of the soliton in the expansion phase exceeds the computational domain.However, since the performed numerical experiments did not demonstrate any noticeable qualitative differences between parametric resonances of solitons with one and two polarizations, this subject was not considered in this work.
However, significant differences were revealed in other processes such as in collisions of two solitons (we recall that the center of each of them moves on a classical tra-  jectory so that they can closely approach in a certain interval of the variable ζ).Collisions of solitons with opposite circular polarizations are significantly different than collisions of identically polarized solitons.In particular, if the left and right solitons with approximately equal numbers N 1 and N 2 approach each other, the first soliton "sees" the second soliton as an unfilled deep potential well for it.Naturally, the first "fluid" tends to "flow" to this second well.Entering the second well, the first soliton transfers a part of its momentum to it.Correspondingly, the second fluid flows to the first well and both wells change their shapes.The final result depends on the relative velocity of longitudinal approach v and on the "impact parameter," i.e., the radius R in terms of strictly helical orbits of solitons in the axisymmetric potential.Thus, even the minimum "set" includes four "input" parameters N 1 , N 2 , v and R. To "trace" the fourdimensional parametric region in sufficient detail, a significant computer time is required.Only collisions with N 1 = N 2 = 5.7 (nonlinearity is not already weak, but the soliton is still stable; the initial values a 2 = b 2 = 0.9 and c 2 = 5.0 are close to equilibrium) have been already simulated more or less systematically.Three values 0.1, 0.2, and 0.3 were taken for the parameter v and the radius R was varied from 0.0 to 1.4 with a step of 0.2.Some characteristic results are presented in Figs.2-4. Figure 2 shows the successive stages of the head-on collision of two solitons with opposite circular polarizations resulting in the formation of two binary light bullets.Figure 3 presents the case R = 1.2 and v = 0.1, where tangential collision led to the merging of two solitons into a more complex rotating nonaxisymmetric structure.Figure 4 (where R = 1.0 and v = 0.1) demonstrates the collision with almost total reflection, where each soliton changed the direction of its motion to opposite, acquiring only a small fraction of the other component.This situation differs from that presented in Fig. 2, where each soliton predominantly continues its motion in the unchanged direction.Some collisions resulted in solitons with approximately equal fractions of both components (e.g., at R = 0.4 and v = 0.1; this case is not shown).
Wave collapse occurred quite often, particularly after secondary collisions near the edge of the computational domain (e.g., at R = 0.2 and v = 0.1; this case is not Several numerical experiments with "heavy" solitons at N 1 = N 2 = 8.0 were also carried out.Such solitons have a low stability margin; therefore, collisions led to collapse in almost all cases except for quite large radii R 1.4.This phenomenon is not discussed here. The approach velocity was low in all reported cases, and the processes of collision were rather long-term and strongly inelastic.The collision of two fast light bullets will be quasi-elastic.Collisions of solitons with the identical circular polarizations were also simulated for comparison.In each of these cases, qualitative differences from the examples reported above were observed, and the development of events strongly depended on the initial phase difference between single-component solitons.For example, the head-on collision such as that presented in Fig. 2 re-sulted in the immediate collapse at zero phase difference, in the repulsion of solitons before close approach at a phase difference of π, and in the flow of the "light fluid" from one soliton to the other at a phase difference of π/2 (these results are not shown).The interaction of solitons with opposite circular polarizations is independent of the phase.

Conclusions
To summarize, the interaction between two polarizations of light has been taken into account for the first time when theoretically considering three-dimensional solitons in graded-index multimode waveguides with the Kerr nonlinearity.The possibility of two-component light bullets has been shown.The numerical simulation of collisions of light bullets with opposite circular polarizations has demonstrated interesting diverse nonlinear dynamics significantly different from the dynamics of collisions of identically polarized solitons.The possible scenarios of the development of events at the scattering of solitons are rich and complex because of a larger number of the internal degrees of freedom that are excited in the two-component system and interact with the translational motion (similar to the one-dimensional case [24][25][26]).Only several first steps have been already made in the study of such binary structures.This field is certainly very promising, in particular, in experiments because single-component solitons in the considered systems have been already implemented in laboratories [5].
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from with the symmetric matrix Ŝ(ζ), where only transverse components are nonzero, the transverse vector F(ζ), and the scalar U 0 (ζ).The functions F(ζ) and U 0 (ζ) include possible deviations of the central line of the waveguide from the axis ζ and variations of the depth of the potential well.

Figure 3 :
Figure 3: Solitons "merged" after the tangential collision at R = 1.2 and v = 0.1.The (x, y) map of the intensity I1(x, y) in the τ = 0 transverse plane at several propagation distances is shown.In this plane, I2(x, y) = I1(−x, −y) due to the symmetry of the initial conditions.

"Figure 4 :
Figure4: Tangential collision at R = 1.0 and v = 0.1 resulting in the formation of two binary solitons.The (τ, x) map of the intensity I1(τ, x) in the y = 0 plane at several propagation distances is shown.In this plane, I2(τ, x) = I1(−τ, −x) due to the symmetry of the initial conditions.Light bullets before and after the collision move along helical lines; for this reason, ζ coordinates at which the centers of solitons are close to the intersection of the y = 0 plane have been chosen.