Elastic collision and molecule formation of spatiotemporal light bullets in a cubic-quintic nonlinear medium S. K. Adhikari∗1 1 Instituto de F́ısica Teórica, UNESP - Universidade Estadual Paulista, 01.140-070 São Paulo, São Paulo, Brazil We consider the statics and dynamics of a stable, mobile three-dimensional (3D) spatiotemporal light bullet in a cubic-quintic nonlinear medium with a focusing cubic nonlinearity above a critical value and any defocusing quintic nonlinearity. The 3D light bullet can propagate with a constant velocity in any direction. Stability of the light bullet under a small perturbation is established numerically. We consider frontal collision between two light bullets with different relative velocities. At large velocities the collision is elastic with the bullets emerge after collision with practically no distortion. At small velocities two bullets coalesce to form a bullet molecule. At a small range of intermediate velocities the localized bullets could form a single entity which expands indefinitely leading to a destruction of the bullets after collision. The present study is based on an analytic Lagrange variational approximation and a full numerical solution of the 3D nonlinear Schrödinger equation. PACS numbers: 05.45.-a, 42.65.Tg, 42.81.Dp I. INTRODUCTION A bright soliton is a self-bound object that travels at a constant velocity in one dimension (1D), due to a cancel- lation of nonlinear attraction and defocusing forces [1, 2]. The 1D soliton in a cubic Kerr medium has been observed in nonlinear optics [1, 2] and in Bose-Einstein conden- sates [3]. Specifically, optical temporal [4] and spatial [5] solitons were observed for a cubic Kerr nonlinearity. However, a three-dimensional (3D) spatiotemporal soli- ton cannot be formed in isolation with a cubic Kerr non- linearity due to collapse [1, 6]. The same is true about a two-dimensional (2D) spatial soliton with a Kerr nonlin- earity. However, the solitons can be stabilized in higher dimensions for a saturable [6, 7] or a modified nonlinear- ity [8], or by a nonlinearity [9] or dispersion [10] manage- ment among other possibilities [11, 12]. A 2D spatiotem- poral optical soliton has been observed [13] in a saturable nonlinearity generated by the cascading of quadratic non- linear processes. The generation of a 2D spatial soliton in an attractive cubic and repulsive quintic medium has been suggested [14] and realized experimentally [15]. The generation of a stable 2D vortex soliton in a cubic-quintic medium has been suggested [16]. Light bullets [6] are localized 3D pulses of electromag- netic energy that can travel through a medium and re- tain their spatiotemporal shape due to a balance between the non-linear self-focusing and spreading effects of the medium in which the pulse beam propagates. Such light bullets are unstable and collapse in a cubic Kerr medium. Light bullets were realized experimentally in arrays of wave guides [17]. There were many theoretical − numer- ical and analytical− studies which established robustness and approximate solitonic nature of the light bullets us- ∗adhikari@ift.unesp.br; URL: http://www.ift.unesp.br/users/adhikari ing the 3D nonlinear Schrödinger (NLS) equation [1] with a modified nonlinearity [7, 8], dissipation [18], and/or dispersion [10]. A saturable nonlinearity leads to stable optical bullets [7]. Nonlinear dissipation in the complex cubic-quintic Ginzburg-Landau equation also stabilizes the bullets [18]. Dispersion management can stabilize light bullets in a medium with cubic nonlinearity [19]. There has been variational study of light bullets in a cubic-quintic medium [20] where a condition of stability was obtained. Another study suggested a way of the sta- bilization of light bullets in a cubic-quintic medium by a periodic variation of diffraction and dispersion [21]. In this paper we study the formation of a 3D spa- tiotemporal light bullet [6, 7] in a cubic-quintic medium for a defocusing quintic nonlinearity and a focusing cu- bic nonlinearity as the ground state of the 3D NLS equa- tion. A cubic-quintic medium is of experimental inter- est also. The study with a polydiacetylene paratoluene sulfonate crystal in the wavelength region near 1600 nm shows that the refractive index versus input intensity cor- relation leads to a cubic-quintic form of nonlinearity in the NLS equation [1, 22]. The cubic-quintic nonlinearity also arises in a low intensity expansion of the saturable nonlinearity used in the pioneering study of light bul- lets [7]. In this study of light bullets in a cubic-quintic medium we find that a stable light bullet can be formed for the cubic focusing nonlinearity above a critical value in the 3D NLS equation for any finite quintic defocusing nonlinearity. The statical properties of the light bullet is studied using a Lagrange variational analysis and a com- plete numerical solution of the 3D NLS equation. The variational and numerical results are found to be in good agreement with each other. The stability of the light bul- let is established numerically under a small perturbation introduced by changing the cubic nonlinearity by a small amount, while the bullet is found to execute sustained breathing oscillation. The light bullet can move freely without deformation along any direction with a constant velocity. We also ar X iv :1 70 1. 03 76 3v 1 [ co nd -m at .q ua nt -g as ] 1 3 Ja n 20 17 http://www.ift.unesp.br/users/adhikari 2 study the frontal collision between two light bullets. Only the collision between two integrable 1D solitons is truly elastic [1]. As the dimensionality of the soliton is in- creased such collision is expected to become inelastic with loss of energy in 2D and 3D. In the present numerical simulation of frontal collision between two light bullets in different parameter domains of nonlinearities and ve- locities three distinct scenarios are found to take place. At sufficiently large velocities the collision is found to be quasi elastic when the two bullets emerge after collision with practically no deformation. At small velocities the collision is inelastic and the bullets form a single bound entity in an excited state and last for ever and execute oscillation. We call this a bullet molecule. In a small domain of intermediate velocities, the bullets coalesce to form a single entity, which expands indefinitely leading to the destruction of the bullets. . We present the 3D NLS equation used in this study in Sec. II. In Sec. III we present the numerical results for stationary profiles of 3D spatiotemporal light bullets. We present numerical tests of stability of the light bul- let under a small perturbation. The quasi-elastic nature of collision of two solitons at large velocities and bullet molecule formation at low velocities are demonstrated by real-time simulation. We end with a summary of our findings in Sec. IV. II. NONLINEAR SCHRÖDINGER EQUATION: VARIATIONAL FORMULATION In nonlinear fibre optics a general 3D NLS equation can usually be written as [1, 23][ i ∂ ∂z + 1 2β0 ∇2 ⊥ + β2 2 ∂2 ∂t2 + γ|A|2 ] A(x, y, t) = 0, (1) ∇2 ⊥ = ∂2 ∂x2 + ∂2 ∂y2 , (2) where β2 is dispersion parameter and can be positive or negative with magnitude of the order of 10−3 ∼ 10−2 ps2/m [23], the nonlinear parameter γ has unit W−1m, the unit of |A|2 is Wm−2, β0 = 2πn0/λ is the prop- agation parameter [1], where n0 is the refractive index and λ is the wave length of the beam. The function A describes the evolution of the beam envelop. For a spa- tiotemporal soliton it is useful to define the characteris- tic lengths for dispersion LDS ≡ τ2/|β2|, and diffraction LDF ≡ 2n0πρ 2/λ where ρ is the radius of the beam, and τ is the time scale of the soliton [11]. For an equi- labrated propagation of the spatiotemporal soliton these two lengths are to be equal − LDS = LDF ≡ LD − yielding ρ2 = LDλ/(2πn0). Now by scaling we define the following dimensionless variables [23] x = x ρ , y = y ρ , t = t τ , z = z LD , φ = Aρ√ P 0 , p = γP0LD ρ2 . (3) The scale P0 is chosen, so that ∫ |φ|2dxdydt = 1. Using the dimensionless variables we obtain the following di- mensionless NLS equations with self-focusing cubic and self-defocusing quintic nonlinearity [1][ i ∂ ∂z + 1 2 ( ∂2 ∂x2 + ∂2 ∂y2 + ∂2 ∂t2 ) + p|φ|2 − q|φ|4 ] φ(r, z) = 0, (4) in scaled units where r ≡ {x, y, t}, p is the cubic and q the quintic nonlinearity. Here z is the propagation distance, x, y denote transverse extensions, and t the time. The plus sign before |φ|2 in Eq. (4) denotes a self-focusing cubic nonlinearity. The quintic nonlinearity of strength q with a negative sign denote self-defocusing. To have an idea of the length and time scales con- cerned, let us consider the case of an infrared beam of wave length λ = 1 µm, and take a nonlinear medium of β2 = 10−2 ps2/m, and a time scale τ = 60 fs. Then the propagation length LD = 36 cm and the beam width ρ ≈ 239 µm. These numbers are quite similar to those in an experiment on spatiotemporal soliton in a planar glass waveguide [24]. In this paper we quote the results in dimensionless units, which can be easily converted to actual experimental units following these guidelines. For an analytic understanding we consider the La- grange variational formulation of the formation of a light bullet. In this spherically symmetric problem, convenient analytic Lagrangian variational approximation of Eq. (4) can be obtained with the following Gaussian ansatz for the wave function [25] φ(r, z) = π−3/4 w3/2(z) exp [ − r2 2w2(z) + iα(z)r2 ] , (5) where r2 = x2 + y2 + t2, w(z) is the width and α(z) is the chirp. The Lagrangian density corresponding to Eq. (4) is given by L(r, z) = i 2 [ φ(r, z) ∂φ∗(r, z) ∂z − φ∗(r, z)∂φ(r, z) ∂z ] + |∇φ(r, z)|2 2 − p 2 |φ(r, z)|4 + q 3 |φ(r, z)|6. (6) Consequently, the effective Lagrangian L ≡ ∫ L(r, z)dr becomes L = 3 2 w2α̇+ 3 4w2 + 3w2α2 − pπ−3/2 4 √ 2w3 + qπ−3 9 √ 3w6 , (7) where the overhead dot denotes the z derivative. The actual physical dimension J/m of this dimensionless La- grangian L can be restored upon multiplication by the factor τP0/LD. The Euler-Lagrange equation for this La- grangian yields the following ordinary differential equa- tion for the width w: ẅ = 1 w3 − p(2π)−3/2 w4 + 4qπ−3 9 √ 3w7 . (8) 3 The energy of the stationary bullet is the Lagrangian (7) with α = 0, e. g., E = 3 4w2 − pπ−3/2 4 √ 2w3 + qπ−3 9 √ 3w6 . (9) The width w of a stationary bullet is obtained by setting the right-hand-side of Eq. (8) to zero: 1 w3 − p(2π)−3/2 w4 + 4qπ−3 9 √ 3w7 = 0, (10) which is the condition for a minimum of energy E − dE/dw = 0, d2E/dw2 > 0 . Without the quintic term (q = 0) the bullet of width w = p/(2π)3/2 is tantamount to an unstable Towne’s soliton [26]. We will demonstrate that for stability a non-zero quintic term (q 6= 0) is nec- essary. For q > 0, Eq. (10) has solution for the cubic nonlinearity p above a critical value pcrit, which is the threshold for the formation of the bullet. III. NUMERICAL RESULTS Unlike in the 1D case, the 3D NLS equation (4) does not have analytic solution and different numerical meth- ods, such as split-step Crank-Nicolson [27] and Fourier spectral [28] methods, are used for its solution. Here we solve it numerically by the split-step Crank-Nicolson method in Cartesian coordinates using a r = {x, y, t} step of 0.025 and a z step of 0.0002 [27]. The number of r discretization points for each components is 256. There are different C and FORTRAN programs for solving the NLS-type equations [27, 29] and one should use the ap- propriate one. We use both imaginary- and real-z prop- agations [27] in the numerical solution of the 3D NLS equation. The imaginary-z propagation is appropriate to find the stationary lowest-energy profile of the bullet. This method replaces z by a new variable ẑ ≡ iz and consequently Eq. (4) becomes completely real and a z- iteration of this equation leads to the lowest-energy state with high accuracy. The real-z propagation involve com- plex variable and hence is more complicated and less ac- curate. However, the real-z propagation yields the prop- agation dynamics of the bullet. In the imaginary-z prop- agation, as the propagation variable z is replaced by the (unphysical) variable ẑ, this method cannot lead to the propagation dynamics of the bullet. In the imaginary-z propagation the initial state was taken as in Eq. (5) with α(z) = 0 and the width w set equal to the variational so- lution obtained by solving Eq. (10). The convergence will be quick if the guess for the width w is close to the final width. All stationary profiles of the bullets are cal- culated by imaginary-z propagation. The dynamics and collision are then studied by real-z propagation using the initial profile obtained in the imaginary-z propagation. The stable bullet corresponds to an energy (E) mini- mum as given by Eq. (10). In Figs. 1(a) and (b) we plot -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (a) q = 30 p = 20 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (a) q = 30 = 30 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (a) q = 30 = 60 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (a) q = 30 = 100 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (a) q = 30 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (b) p = 60 q = 10 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (b) p = 60 = 30 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (b) p = 60 = 60 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (b) p = 60 = 100 -40 -30 -20 -10 0 10 0 0.5 1 1.5 2 E (w ) w (b) p = 60 0 10 20 30 0 20 40 60 80 100 p c ri t q (c) -4 -3 -2 -1 0 1 2 0 0.5 1 1.5 2 E (w ) w (d) p = 10 q = 0.6 -4 -3 -2 -1 0 1 2 0 0.5 1 1.5 2 E (w ) w (d) p = 10 q = 0.6 = 0.8 -4 -3 -2 -1 0 1 2 0 0.5 1 1.5 2 E (w ) w (d) p = 10 q = 0.6 = 0.8 = 1.3 -4 -3 -2 -1 0 1 2 0 0.5 1 1.5 2 E (w ) w (d) p = 10 q = 0.6 = 0.8 = 1.3 = 2 -4 -3 -2 -1 0 1 2 0 0.5 1 1.5 2 E (w ) w (d) p = 10 q = 0.6 = 0.8 = 1.3 = 2 = 5 -4 -3 -2 -1 0 1 2 0 0.5 1 1.5 2 E (w ) w (d) p = 10 FIG. 1: (Color online) (a) Variational energy versus width (E − w) plot for different cubic nonlinearities p(= 20, 30, 60, 100) and quintic nonlinearity q = 30. (b)The same for different quintic nonlinearities q (= 10, 30, 60, 100) and cu- bic nonlinearity p = 60. (c) Variational critical cubic nonlin- earity pcrit for light bullet formation, obtained from Eq. (10), for different values of quintic nonlinearity q. (d) Variational energy versus width (E −w) plot for different quintic nonlin- earities q(= 0.6, 0.8, 1.3, 2, 5) and cubic nonlinearity p = 10. E versus w of Eq. (10) for different cubic (p) and quin- tic (q) nonlinearities. The energy minima of these plots correspond to a stable bullet of negative energy. From Fig. 1(a) we find that for q = 30 such an energy minima exists for p > 20. An accurate value of this limit can be obtained from Eq. (10): for q = 30 this equation has solution for p ≥ pcrit = 19.6. Hence the NLS equation 4 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 r rm s p (a) var (q = 10) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 r rm s p (a) var (q = 10) num (= 10) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 r rm s p (a) var (q = 10) num (= 10) var (q = 30) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 r rm s p (a) var (q = 10) num (= 10) var (q = 30) num (= 30) 0 20 40 60 80 100 0 20 40 60 80 100 |E | p (b) var (q = 10) 0 20 40 60 80 100 0 20 40 60 80 100 |E | p (b) var (q = 10) num (= 10) 0 20 40 60 80 100 0 20 40 60 80 100 |E | p (b) var (q = 10) num (= 10) var (q = 30) 0 20 40 60 80 100 0 20 40 60 80 100 |E | p (b) var (q = 10) num (= 10) var (q = 30) num (= 30) FIG. 2: (Color online) Variational (var) and numerical (num) (a) rms radius and (b) energy |E| versus cubic nonlinearity p of a light bullet for two different quintic nonlinearities q (= 10, 30). (4) can have a stable light bullet solution for cubic non- linearity p greater than a critical value pcrit. For p < pcrit the system is much too repulsive and is not bound and escapes to infinity. However, this critical value pcrit of p is a function of the quintic nonlinearity q. The pcrit−q cor- relation can be found from an attempt to solve Eq. (10) numerically. The pcrit − q correlation obtained in this fashion is plotted in Fig. 1(c). However, in addition to these stable bullets corresponding to a global minimum of energy with negetive energy values, one can also have metastable bullets corresponding to a local minimum of energy at positive energies. Such a situation is illustrated in Fig. 1(d) where we plot the variational E − w curves for p = 10 and different q values. The bullet with p = 10 and q = 1.3 has a local minimum at a positive energy and is metastable in nature. In the following we will only consider the stable light bullets with negative energy. Next we compare in Fig. 2(a) the numerical and varia- tional root-mean-square (rms) radius rrms of a light bul- let versus cubic nonlinearity p for different quintic non- linearity q = 10, 30. The variational result is given by: rrms = √ 3/2w, where w is the equilibrium variational width. In Fig. 2(b) we show the numerical and varia- tional energy |E| of a light bullet versus p for different q. The numerical energy is calculated using Eq. (9) with the numerically obtained φ(r, z). The energy of the light 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x p = 35, q = 2 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x p = 35, q = 2 = 30, = 5 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x p = 35, q = 2 = 30, = 5 = 50, = 15 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x p = 35, q = 2 = 30, = 5 = 50, = 15 = 30, = 15 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x p = 35, q = 2 = 30, = 5 = 50, = 15 = 30, = 15 = 20, = 15 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x p = 35, q = 2 = 30, = 5 = 50, = 15 = 30, = 15 = 20, = 15 = 25, = 35 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x 0 1 2 3 -2 -1 0 1 2 ρ 1 D (x ) x FIG. 3: (Color online) Numerical (line) and variational (chain of symbols) reduced 1D density ρ1D(x) for different cubic non- linearity p and quintic nonlinearity q. bullet is negative in all cases and its absolute value is plot- ted. For small p, the agreement between numerical and variational results is better. For large p, the agreement between the two is qualitative. For large nonlinearity p in the NLS equation, the profile of the bullet deviates more from the Gaussian variational ansatz − thus making the variational results more approximate. To study the density distribution of the light bullets we calculate the reduced 1D density defined by ρ1D(x) = ∫ dtdy|φ(r)|2. (11) In Fig. 3 we plot this reduced 1D density as obtained from variational and numerical calculations for different cubic nonlinearity p and quintic nonlinearity q. For a fixed defocusing nonlinearity q (= 15), the light bullet is more compact with the increase of focusing nonlinearity p resulting in more attraction. For a fixed focusing non- linearity p (= 20), the light bullet is more compact with the decrease of defocusing nonlinearity q resulting in less repulsion. Now we present a numerical test of stability of a sta- ble bullet under a small perturbation. For this purpose we consider the bullet shown in Fig. 3 with p = 20 and q = 15 as calculated by imaginary-z propagation. Using the imaginary-z profile as the initial state we perform nu- merical simulation by real-z propagation under a small perturbation introduced at z = 0 by changing p from 20 to 20.5. This sudden perturbation in the cubic non- linearity increases the attraction in the system and the light bullet starts a rapid breathing oscillation. In Fig. 4 we show the steady oscillation in the rms x size xrms versus propagation distance z during real-z propagation. The steady continued oscillation of the bullet over a long distance of propagation establishes the stability of the bullet. The real-z simulation was performed in full 3D space without assuming spherical symmetry to guaranty the stability in full 3D Cartesian space. The collision between two analytic 1D solitons is truly elastic [1, 3] and such solitons pass through each other 5 0.4 0.41 0.42 0.43 0 10 20 30 40 x rm s z FIG. 4: (Color online) Steady oscillation of the rms size xrms during real-z propagation of a light bullet with p = 20 and q = 15, when at z = 0 the cubic nonlinearity p is suddenly changed from 20 to 20.5. (a) 0 0.06 0.12 0.18 z -5 -4 -3 -2 -1 0 1 2 3 4 5 x 0 1 2 3 4 ρ1D(x,z) (c) 0 0.06 0.12 0.18 z -5 -4 -3 -2 -1 0 1 2 3 4 5 x 0 1 2 3 4 ρ1D(x,z) FIG. 5: (Color online) (a) The 1D density ρ1D(x, z) and (b) its contour plot during the collision of two light bullets of Fig. 3 with p = 20, q = 15 initially placed at x = ±3 , upon real- z propagation. The initial wave functions are multiplied by exp(±i40x) which sets them in motion with velocity of about 33 each. The same for two light bullets with p = 25, q = 35 of Fig. 3 are shown in (c) and (d). without deformation at any incident velocities. The col- lision between two 3D spatiotemporal light bullets is ex- pected to be inelastic in general with loss of kinetic en- ergy resulting in the deformation of the bullets. Such a collision can at best be quasi-elastic. To test the solitonic nature of the present light bullets, we study the frontal head-on collision of two bullets. The imaginary-z profiles of the light bullets shown in Fig. 3 with (i) p = 20, q = 15 and (ii) p = 25, q = 35 are used as the initial function in the real-z simulation of collision, with two identical bullets placed at x = ±3 initially for z = 0. To set the light bullets in motion along the x axis in opposite FIG. 6: (Color online) 3D isodensity profile of the (a) initial (at z = 0) and (b) final (at z = 0.18) light bullets each with p = 20, q = 15 and (c) initial (at z = 0) and (d) final (at z = 0.18) light bullets each with p = 25, q = 35 undergoing elastic collision illustrated in Fig. 5(a) and (b). density on contour 0.01. directions the respective wave functions are multiplied by exp(±i40x). To illustrate the dynamics upon real- z simulation, we plot the time evolution of 1D density ρ1D(x, z) ≡ ∫ dt ∫ dy|φ(r, z)|2 in Fig. 5(a) for the colli- sion of two bullets with p = 20, q = 15. The correspond- ing contour plot is presented in Fig. 5(b). The same for the collision of two light bullets with p = 25, q = 35 is presented in Figs. 5(c) and (d). The dimensionless velocity of a light bullet is about 33 and the deviation from elastic collision is found to be small. Considering the three-dimensional nature of collision, the distortion in the bullet profile is found to be negligible. To visu- alize the amount of inelasticity in the collision displayed in Figs. 5(a) and (b), we display in Figs. 6(a) and (b) the 3D isodensity profile of the light bullet before (z = 0) and after (z = 0.18) the collision shown in Figs. 5(a) and (b). The same for the collision shown in Figs. 5(c) and (d) is illustrated in Figs. 6(c) and (d). In general, the inelasticity in both collisions is small. To study the inelastic collision we consider two com- pact bullets with p = 35, q = 2 and place them at x = ±1 and set them in motion with dimensionless velocity of about 2 in opposite directions. This is achieved by mul- tiplying the respective wave functions by exp(±i2x) and perform real-z simulation. The dynamics is illustrated by a plot of the time evolution of 1D density ρ1D(x, z) in Fig. 7 (a) and the corresponding contour plot is shown in Fig. 7 (b). The same for the collision of two light bul- lets with p = 30, q = 15 with an initial velocity of about 0.5 are illustrated in Figs. 7(c) and (d), respectively. In 6 (a) 0 1 z -2 -1 0 1 2 x 0 5 10 ρ1D(x,z) (c) 0 5 10 z -4 -2 0 2 4 x 0 2 4 ρ1D(x,z) FIG. 7: (Color online) (a) The 1D density ρ1D(x, z) and (b) its contour plot during the collision of two light bullets of Fig. 3 with p = 35, q = 2 initially placed at x = ±1, upon real- z propagation. The initial wave functions are multiplied by exp(±i2x) to set them in motion with an initial dimensionless velocity of about 2. The same for the collision of two light bullets with p = 30, q = 15 of Fig. 3 for an initial velocity of about 0.5 are shown in (c) and (d), respectively. both cases the two bullets come close to each other at x = 0 coalesce to form a bullet molecule and never sepa- rate again. The combined bound system remain at rest at x = 0 continuing small breathing oscillation because of a small amount of liberated kinetic energy which creates the bullet molecule in an excited state. The observation of oscillating bullet molecules has been reported some time ago in dissipative systems [30]. Hence at sufficiently small incident velocities the colli- sion of two light bullets lead to the formation of a bullet molecule and at large velocities one has the quasi-elastic (a) 0 0.4 0.8 1.2 1.6 z -9 -6 -3 0 3 6 9 x 0 1 2 3 4 ρ1D(x,z) FIG. 8: (Color online) (a) The 1D density ρ1D(x, z) and (b) its contour plot during the collision of two light bullets of Fig. 3 with p = 30, q = 15 initially placed at x = ±4.5, upon real-z propagation. The initial wave functions are multiplied by exp(±i6x) to set them in motion with an initial velocity of about 6. collision of two light bullets. At intermediate velocities a new phenomenon can take place. As the initial velocity is slowly increased from the domain of molecule formation, after collision a bullet molecule is formed in a highly ex- cited state with a large amount of energy. In that case, because of the excess energy the bullet molecule expands and the localized bullets are destroyed. This is illustrated by a plot of the time evolution of 1D density ρ1D(x, z) in Fig. 8 (a) for the case of collision of two bullets with p = 30, q = 15 at an initial velocity of about 6 and the corresponding contour plot is shown in Fig. 8 (b). IV. SUMMARY Summarizing, we demonstrated the creation of a stable 3D spatiotemporal light bullet with cubic-quintic nonlin- earity employing the Lagrange variational and full 3D numerical solution of the NLS equation. The statical properties of the bullet are studied by a variational ap- proximation and a numerical imaginary-z solution of the 3D NLS equation. The cubic nonlinearity is taken as focusing Kerr type above a critical value, whereas the quintic nonlinearity is defocusing. The dynamical prop- erties are studied by a real-z solution of the NLS equa- tion. In the 3D spatiotemporal case, the optical bullet can move with a constant velocity. At large velocities, the collision between the two spatiotemporal light bul- lets is quasi elastic with no visible deformation of the final bullets. At small velocities, the collision is inelastic with the formation of a bullet molecule after collision. At medium velocities the bullets can be destroyed after collision. Acknowledgments I thank F. K. Abdullaev for very valuable comments. Interesting discussion with Boris A. Malomed is grate- fully acknowledged. We thank the Fundação de Amparo à Pesquisa do Estado de São Paulo (Brazil) (Project: 2012/00451-0 and the Conselho Nacional de Desen- volvimento Cient́ıfico e Tecnológico (Brazil) (Project: 303280/2014-0) for support. 7 [1] Y. S. Kivshar and G. Agrawal, Optical Solitons: From Fibers to Photonic Crystals, (Academic Press, San Diego, 2003). [2] Y. S. Kivshar and B. A. Malomed, Rev. Mod. Phys. 61, 763 (1989); V. S. Bagnato, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and D. Mihalache, Rom. Rep. Phys. 67, 5 (2015); D. Mihalache, Rom. J. Phys. 59, 295 (2014). [3] V. M. Perez-Garcia, H. Michinel, and H. Herrero, Phys. Rev. A 57, 3837 (1998); F. K. Abdullaev, A. Gammal, A. M. Kamchatnov, and L. Tomio, Int. J. Mod. Phys. B 19, 3415 (2005); K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet, Nature 417, 150 (2002); L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cu- bizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 256, 1290 (2002). [4] P. Di Trapani, D. Caironi, G. Valiulis, A. Dubietis, R. Danielius, and A. Piskarskas, Phys. Rev. Lett. 81, 570 (1998). [5] J. U. Kang, G. I. Stegeman, and J. S. Aitchison, Opt. Lett. 21, 189 (1996); J. U. Kang, J. S. Aitchison, G. I. Stegeman, and N. N. Akhmediev, Opt. Quantum Elec- tron. 30, 649 (1998); J. U. Kang, G. I. Stegeman, J. S. Aitchison, N. Akhmediev, Phys. Rev. Lett. 76, 3699 (1996). [6] Y. Silberberg, Opt. Lett. 15, 1282 (1990). [7] N. Akhmediev and J. M. Soto-Crespo, Phys. Rev. A 47, 1358 (1993); D. V. Skryabin and W. J. Firth, Opt. Com- mun. 148, 79 (1998); D. E. Edmundson and R. H. Enns, Opt. Lett. 17, 586 (1992); D. E. Edmundson and R. H. Enns, Opt. Lett. 18, 1609 (1993); G. Fibich and B. Ilan, Opt. Lett. 29, 887 (2004); L. Torner and Y. V. Kartashov, Opt. Lett. 34, 1129 (2009); D. E. Edmund- son and R. H. Enns, Phys. Rev. A 51, 2491 (1995); A. A. Kanashov and A. M. Rubenchik, Physica D 4, 122 (1981); D. Mihalache, D. Mazilu, L.-C. Crasovan, I. Tow- ers, A. V. Buryak, B. A. Malomed, L. Torner, J. P. Tor- res, and F. Lederer, Phys. Rev. Lett. 88, 073902 (2002). [8] R. McLeod, K. Wagner, and S. Blair, Phys. Rev. A 52, 3254 (1995); D. Mihalache, D. Mazilu, L.-C. Crasovan, L. Torner, B. A. Malomed, and F. Lederer, Phys. Rev. E 62, 7340 (2000); O. Bang, W. Krolikowski, J. Wyller, and J. J. Rasmussen, Phys. Rev. E 66, 046619 (2002). [9] S. K. Adhikari, Phys. Rev. A 69, 063613 (2004), Phys. Rev. E 70, 036608 (2004); H. Saito and M. Ueda, Phys. Rev. Lett. 90, 040403 (2003); F. K. Abdullaev, J. G. Caputo, R. A. Kraenkel, and B. A. Malomed, Phys. Rev. A 67, 013605 (2003). [10] S. K. Adhikari, Phys. Rev. E 71, 016611 (2005). [11] B. A. Malomed, D. Mihalache, F. Wise, and L. Torner, J. Opt. B 7, R53 (2005). [12] B. A. Malomed, F. Lederer, D. Mazilu, and D. Mihalache, Phys. Lett. A 361, 336 (2007). [13] X. Liu, L. J. Qian, and F. W. Wise, Phys. Rev. Lett. 82, 4631 (1999). [14] M. L. Quiroga-Teixeiro, A. Berntson, and H. Michinel, J. Opt. Soc. Am. B16, 1697 (1999). [15] E. L. Falcão-Filho, C. B. de Araujo, G. Boudebs, H. Leblond, and V. Skarka, Phys. Rev. Lett. 110, 013901 (2013). [16] V. I. Berezhiani, V. Skarka, and N. B. Aleksić, Phys. Rev. E 64, 057601 (2001). [17] S. Minardi, F. Eilenberger, Y. V. Kartashov, A. Szameit, U. Röpke, J. Kobelke, K. Schuster, H. Bartelt, S. Nolte, L. Torner, F. Lederer, A. Tünnermann, and T. Pertsch, Phys. Rev. Lett. 105, 263901 (2010). [18] B. Liu and X.-D. He, Opt. Express 19, 20009 (2010); P. Grelu, J. Soto-Crespo, and N. Akhmediev, Opt. Express, 13, 9352 (2005); N. B. Aleksić, V. Skarka, D. V. Timo- tijević, and D. Gauthier, Phys. Rev. A 75, 061802(R) (2007);A. Kamagate, Ph. Grelu, P. Tchofo-Dinda, J. M. Soto-Crespo, and N. Akhmediev, Phys. Rev. E 79, 026609 (2009); S. Chen, Phys. Rev. A 86, 033829 (2012); D. Mihalache, D. Mazilu, F. Lederer, Y. V. Kartashov, L. -C. Crasovan, L. Torner, and B. A. Malomed, Phys. Rev. Lett. 97, 073904 (2006) J. M. Soto-Crespo, Ph. Grelu, and N. Akhmediev, Opt. Express 14, 4013 (2006); N. N. Rozanov, J. Opt. Technol. 76, 187 (2009); D. Mihalache, D. Mazilu, F. Lederer, H. Leblond, and B. A. Malomed, Eur. Phys. J. Special Topics 173, 245 (2009). [19] M. Matuszewski, E. Infeld, B. A. Malomed, and M. Trip- penbach, Opt. Commun. 259, 49 (2006). [20] Z. Jovanoski, J. Mod. Opt. 48, 865 (2001). [21] P. A. Subha, K. K. Abdullah, and V. C. Kuriakose, J. Nonlinear Opt. Phys. Mat. 19, 459 (2010). [22] B. Lawrence, W. E. Torruellas, M. Cha, M. L. Sund- heimer, G. I. Stegeman, J. Meth, S. Etemad, and G. Baker, Phys. Rev. Lett. 73, 597 (1994). [23] G. P. Agrawal, J. Opt. Soc. Am. B 28, A1 (2011). [24] H. S. Eisenberg, R. Morandotti, Y. Silberberg, S. Bar- Ad, D. Ross, and J. S. Aitchison, Phys. Rev. Lett. 87, 043902 (2001). [25] V. M. Perez-Garcia, H. Michinel, J. I. Cirac, M. Lewen- stein, and P. Zoller, Phys. Rev. Lett. 77, 5320 (1996). [26] R. Y. Chiao, E. Garmire, and C. H. Townes, Phys. Rev. Lett. 13, 479 (1964). [27] P. Muruganandam and S. K. Adhikari, Comput. Phys. Commun. 180, 1888 (2009); J. Phys. B 36, 2501 (2003); D. Vudragović, I. Vidanović, A. Balaž, P. Muruganan- dam, and S. K. Adhikari, Comput. Phys. Commun. 183, 2021 (2012); L. E. Young-S., D. Vudragović, P. Muru- ganandam, S. K. Adhikari, and A. Balaž, Comput. Phys. Commun. 204, 209 (2016). [28] P. Muruganandam and S. K. Adhikari, J. Phys. B 36, 2501 (2003). [29] B. Satarić, V. Slavnić, A. Belić, Antun Balaž, P. Muru- ganandam, and S. K. Adhikari, Comput. Phys. Commun. 200, 411 (2016); V. Loncar, A. Balaž, A. Bogojević, S. Skrbić, P. Muruganandam, and S. K. Adhikari, Comput. Phys. Commun. 200, 406 (2016); R. Kishor Kumar, L. E. Young-S., D. Vudragović, Antun Balaž, P. Muruganan- dam, and S. K. Adhikari, Comput. Phys. Commun. 195, 117 (2015). [30] J. M. Soto-Crespo, N. Akhmediev, and Ph. Grelu, Phys. Rev. E 74, 046612 (2006); N. Akhmediev, J. M. Soto- Crespo, and Ph. Grelu, Chaos 17, 037112 (2007). I Introduction II Nonlinear Schrödinger equation: Variational formulation III Numerical Results IV Summary Acknowledgments References