Trapping of two-component light bullets in a gradient waveguide with normal group dispersion

In this paper we consider the process of the second harmonic generation in a gradient waveguide, taking into account diffraction and relatively weak temporal dispersion. Using the slowly varying envelope approximation and neglecting the dispersion of the nonlinear part of the response of the medium we obtain the system of parabolic equations for the envelopes of both harmonics. We also derive integrals of motion of this system. To solve it numerically we construct a nonlinear finite-difference scheme based on the Crank-Nicolson method preserving the integrals. Primarily, we focus our investigations on the processes of a two-component light bullets generation. We demonstrate that the generation of a coupled pair is possible in a planar waveguide even at normal group velocity dispersion.


Introduction
Propagation of waves in a homogeneous boundless medium is an idealization which is rarely found in nature. In most cases, the properties of the medium vary in space and this leads to a significant change in the nature of wave propagation.
Sometimes it is reasonable to add inhomogeneity to reach a desirable scenario of wave propagation. For instance, it is widely known that multi-dimensional light bullets are highly unstable in a homogeneous medium with Kerr nonlinearity. However, if we use an inhomogeneous Kerr medium, where the linear part of the refraction index depends on transverse coordinates, the wave collapse can be prevented [1][2][3][4]. The interplay between dispersion, diffraction, index inhomogeneity and Kerr nonlinearity was considered in detail in these papers. It was shown that, it is possible to form light bullets in a gradient nonlinear medium in the regime of self-defocusing and normal dispersion [3].
Light bullets at quadratic nonlinearity are more attractive than at cubic one [4]: they have a lower excitation threshold, they are more stable. Multi-component optical solitons at quadratic nonlinearity were firstly predicted in 1974 [5] when it was demonstrated that one-dimensional spatial or temporal pulse spreading could be compensated by linear self-compression. To date, the theory of multi-dimensional multi-component quadratic solitons has been essentially developed. PLOS  Spatio-temporal solitons, called also light bullets, were carefully studied either theoretically or experimentally by numerous researchers [4,6,7]. Recently, in [8,9] we presented a detailed theory of breathing light bullets in a homogeneous medium with quadratic nonlinearity at anomalous dispersion. Using the averaged Lagrangian method, we derived approximate analytical solutions in the form of two-component planar spatio-temporal solitons. Besides that, an appropriate physical model was suggested for the THz range and anomalous dispersion at both frequencies.
It is known that stable multi-component optical solitons at quadratic nonlinearity are usually observed under anomalous dispersion. Including normal dispersion in consideration, one may enlarge the range of frequencies of wave localization. However, the problem of simultaneous compensation of either nonlinear or linear effects is a challenging one. By linear effects we mean diffraction and normal dispersion stretching a wave packet. As it has been already mentioned above, waveguides demonstrate a remarkable ability to support multi-dimensional soliton structures in the media with Kerr nonlinearity [1][2][3][4]. Their geometry may be chosen in such a way that it may play either focusing or defocusing role and compete with stretching tendencies. The idea to apply a similar approach to the medium with quadratic nonlinearity seems to be a promising one.
The present study of light bullets in a planar waveguide at quadratic nonlinearity is undertaken as a logical continuation of our previous investigations for homogeneous media [8,9]. A crucial role of the competition between quadratic nonlinearity, dispersion, diffraction and inhomogeneity was discussed shortly in [10,11].
In this paper we carefully discuss physical features of the second harmonic generation process in a gradient waveguide, taking into account diffraction and relatively weak temporal dispersion. Primarily, we focus on the process of the second harmonic generation (SHG) and the birth of a coupled pair (two-component light bullet). Using the slowly varying envelope approximation and neglecting dispersion of the nonlinear part of the medium response, we develop a governing system of equations for SHG in a waveguide with transverse inhomogeneity. We prove this system possesses two motion integrals. Thus, a numerical algorithm used for simulation must preserve difference analogs of these integrals [12,13]. Besides that, we are dealing with a multidimensional problem and, therefore, it is especially important to use algorithms that save time. Usually the splitting technique is applied to this purpose (see [12,14,15] and references therein). But such algorithms are not appropriate if it is necessary to preserve motion integrals of the class of problems to which the considered problem belongs [16]. Two other general approaches to the construction of numerical methods for the Schrӧdinger equations exist, namely, Fast Fourier Transform (FFT) (see also [17]) and multi-grid (MG) techniques [15,17]. Splitting technique, FFT and MG are compared in [15]. Efficiency of MG and FFT are found close to each other. FFT methods were successfully applied to the problems of nonlinear optics [6,18,19]. But even FFT solvers are not sufficiently efficient in the case of 2D +1 or 3D+1 problems. Instead of them, multi-step iterative procedures could be applied for the implementation of multi-dimensional conservative difference scheme. In [20,21] a two-step iterative algorithm is proposed for a model describing the process of femtosecond optical pulse propagation in semiconductors. As the Crank-Nicolson method is used in this approach, it is conservative and is promising in efficiency which is comparable with that of the splitting technique. Thus, for the problem under consideration in our study we construct a conservative nonlinear finite-difference scheme and develop a multi-step iterative algorithm for its implementation.
The paper is organized as follows. In Sec. 2, using the quasi-optical approach, we introduce the governing equations' system and consider in detail the way which we utilize for the description of transversal inhomogeneity. In Sec. 3 we describe the numerical algorithm used to perform direct numerical simulation and discuss its properties. Results of direct numerical simulation are discussed in Sec. 4. Sec. 5 contains the conclusions.

Quasi-optical approach
The way to derive the equation describing the interplay between dispersion, diffraction, inhomogeneity, and nonlinearity was described in detail in [3] for Kerr nonlinearity. There the authors started from Maxwell's equations supplemented by the dependence of the refractive index on carrier frequency and transversal coordinates. Then a standard procedure was fulfilled, the paraxial and the slowly varying envelope approximations were applied. It is important to underline that the equation for a slowly varying envelope obtained in that paper is similar to the standard multidimensional nonlinear Schrodinger equation (NLSE) except the only one additional term due to nonlinear medium inhomogeneity.
In the present work for the consideration of the SHG in a quadratically nonlinear waveguide we applied about the same procedure.
To describe waveguide geometry we represent the linear frequency susceptibility χ ω (r ? ) in the form Here r ? is the transverse radius-vector perpendicular to the central axis of the waveguide, w ð0Þ o is the linear susceptibility of the medium at the center of the waveguide cross section; the dimensionless function f ω (r ? ) characterizes the transverse inhomogeneity of the susceptibility and satisfies the condition f ω (0) = 0. Similar to [3], the paraxial approximation is also applied. We assume that both harmonics propagate along the z-axis, and at that, we suppose the linear group velocities v ðoÞ g and v ð2oÞ g (envelope carrier frequencies ω and 2ω correspondingly) at the center of the waveguide (r ? = 0) are related by the following inequality jv ð2oÞ Provided inhomogeneity, nonlinearity, dispersion, and diffraction are weak, we use the slowly varying envelope approximation [22]. The resulting system of equations for the envelopes of the fundamental F 1 and second F 2 harmonics looks as follows: f ω (r ? ) and f 2ω (r ? ) are the dimensionless functions satisfying the condition f ω (0) = f 2ω (0) = 0, w ð2Þ ðo; oÞ. χ (2) (2ω,−ω) and χ (2) (ω,ω) are the second order nonlinear optical susceptibilities at the waveguide center.
The first and second terms in the right-hand sides of Eq (1.1) describe the effect of the waveguide (transverse inhomogeneity) on the phase and group velocities of the harmonics, respectively. Spatial inhomogeneity of the linear refractive indices for both harmonics in (1.1) is taken into account in the same way as it was done in [3,4] when studying the waveguide propagation mode of the soliton of the NLSE.
Putting in (1.1) g 1,2 = 0, we come to the well-known system of equations for the pulsed mode of SHG in a homogeneous medium [23,24].
Dependences of the refractive indexes n ω,2ω on the transverse coordinates r ? of the waveguide are expressed through the functions f ω,2ω (r ? ) as follows: If we deal with a focusing waveguide, then the functions f ω (r ? ) and f 2ω (r ? ) decrease from the center to the periphery. In the opposite case, the waveguide is defocusing. Below we consider planar waveguides, i.e. Δ ? = @ 2 /@x 2 and f ω,2ω = f ω,2ω (x), the profiles of the waveguide functions are conveniently chosen, for example, in the form of parabolic profile with saturation: Here x is the distance from the centre of the waveguide to a current point in transversal direction, ε w = 1 for the defocusing waveguide and ε w = −1 for the focusing waveguide.
In addition, we can consider profiles in the forms Note that in the case of focusing waveguides (ε w = −1) the transverse susceptibility profiles (see (1.3)) for the dependences (1.4), (1.5) and (1.6) have the Lorentzian, exponential and Gaussian modes, respectively: ð1:7Þ ð1:8Þ ð1:9Þ In all these cases, the refractive indexes n ω,2ω (x) decrease from the maximum value at the center of the waveguide (at x = 0) to unity at its periphery (at x!1).

Dimensionless equations and integrals of motion
We still consider a planar (Δ ? = @ 2 /@x 2 ) waveguide and use dimensionless parameters related to the physical parameters in the following way: Here A in is the input peak amplitude of the fundamental harmonic, R in and τ in are initial pulse spatial and temporal widths, respectively. We introduce also the following propagation and waveguide characteristics: a 2ω are the characteristic lengths of waveguide transversal inhomogeneity. Then, we suppose the group velocity matching conditions are satisfied (v ðoÞ . Finally, we get the following system of the dimensionless equations which we use as a base for our numerical simulation: Boundary and initial conditions are as follows. In (2.1)-(2.2) L z is the dimensionless length of the nonlinear medium, L τ is the dimensionless time interval during which laser pulse interaction with the medium is analyzed. L x is the dimensionless length of the transversal domain. We choose finite lengths of the transversal coordinate and time with zero conditions at the boundaries of these coordinates from the following considerations. Since we deal with a finite pulse, it is naturally to take a sufficiently long time interval whose boundaries are not influenced by the pulse. Besides, the transversal size of the studied bullet is also finite due to narrow laser radiation. Thus, we can choose a wide transversal domain with boundaries which are not affected by radiation. Obviously, this choice assumes the need to monitor the fulfillment of zero boundary conditions, and, if necessary, a widening of the computational domain. It may cause specific computational difficulties, which nevertheless, can be eliminated by imposing artificial boundary conditions (see [25] and references therein). In (2.1) p 1,2 are dimensionless functions describing waveguide inhomogeneity and corresponding to (1.4)-(1.6). Below we discuss simulation with ; ð2:3Þ ; ð2:4Þ ð2:5Þ Here � x is the distance from the center of the waveguide to the current point. In general, a ω and a 2ω are not equal to each other ( 3) describes a waveguide with parabolic profile, (2.4)-(2.5) correspond to profiles with Lorentzian (2.4), tangential (2.5), and exponential (2.6) saturation. Below we will omit bars in the notations of dimensionless variables.
The system (2.1)-(2.2) possesses the motion integrals In (2.8) φ 1,2 are the wave phases. To derive (2.7) we multiply the equations of (2.1) by gA � 1 ; A � 2 correspondingly, then integrate both parts of these transformed equations with respect to x,τ, sum them up and take real parts of all components. Finally we get the equality (2.7) expressing a law of energy conservation. At the next step we proceed, multiplying the equations of (2.1) by 2g @A � 1 @z ; @A � 2 @z correspondingly. Then we repeat the successive steps used when receiving (2.7), but at the last step we take imaginary parts of all components. So, we obtain (2.8) concerning the evolution of the phases of harmonics.

Approximation of the equations
We introduce uniform grids in the domain Γ and in the z domain: ð3:1Þ A numerical approximation to the exact solution of the problem (2.1)-(2.2) A l;j;k 1;2 ¼ A 1;2 ðz l ; x j ; t k Þ we consider on the grid ω z ×ω Γ and denote it by C l;j;k 1;2 ¼ C 1;2 ðz l ; x j ; t k Þ. To approximate first-and second-order derivatives with respect to x and τ we use the standard expressions: C l;j;k 1;2 � tt ¼ Initial and boundary conditions are approximated exactly.
This scheme is known to be of the second order of approximation with respect to all coordinates [12][13][14]. It is easily generalized to the case of phase mismatch (Δk6 ¼0).

Two-step iterative process
Nonlinear scheme (3.2)-(3.3) can be resolved with the help of an iteration process [12]. But the computational complexity of the direct matrix inversion after linearization makes such an approach practically useless [12,15]. FFT technique [15,19] is not straightforward in the case under consideration due to the x-dependent coefficients p 1,2 in (3.2). Therefore, we develop the approach proposed in [20,21] and write down the following two-step iteration process for the implementation of (3.2)-(3.3).
At the first step we seek for the iteration (s+1) of the difference functions C sþ1 1;2 l þ 1; j; k : At the second step we complete the procedure at the current iteration cycle finding C sþ2 1;2 l þ 1; j; k : We see that at each step of the iteration process we deal with one-dimensional problem. Thus

. Provided that h z �C(h x h τ ) 2 the unique solution of the difference scheme (3.2)-(3.3) exists and the two-step iteration process (3.4)-(3.7) converges to it as a geometric progression with the denominator q�h z /(h x h τ ) 2 .
Proof. The proof is made with the help of the contraction mapping theorem. We have to demonstrate that (3.4)-(3.7) is a contraction, i.e that all iterations are uniformly limited and that in a certain difference norm k�k ( we use at each step of the iteration process the Cauchy-Bunyakovsky-Schwarz inequality, estimate the inverse norm of the difference operators F � xx and F � tt [12]. This procedure allows us firstly, to show that the iterations are uniformly limited if h z �C(h x h τ ) 2 . Then, considering consecutively the iteration differences

Conservativeness Theorem 3.2. The scheme (3.1)-(3.2) is conservative: it preserves the following difference analogues of the integrals (2.7)-(2.8) in the norm
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffiffi jC lþ1 ð3:10Þ Proof. We follow the procedure described in [13]. Multiplying both parts of the first equation of (3.1) by ðgC l;j;k

Initial pulse and absorbing boundary conditions
In this Section we discuss the outcome of a direct numerical simulation of the system (2.1) on the base of the finite-difference scheme (3.2). In computations we launch the initial pulse at both frequencies or at the fundamental frequency only Implementing calculations we usually choose L x and L τ so that the amplitudes of both harmonics decay to zero at the boundaries of this domain. In particular cases, when L x and L τ are too large, the efficiency of the method can be improved by using absorbing boundary conditions along the coordinates x and τ. To this end we embed an artificial absorption in the system (2.1) defining the function σ(x,τ) = σ min (σ x (x)+σ τ (τ)), where  Analyzing Fig 1(A), which represents linear diffraction of the pulse, one can see that changing the value of σ min at fixed values of x ab , L x , and ε xτ , we observe the best coincidence with the reference curve obtained with the larger computational transversal domain (black dotted line) at σ min = 0.1 (green line). Fig 1(B) demonstrates the choice of optimal ε xτ = 5 in a similar way. Here we fix x ab , L x , and σ min and vary ε xτ . Parameter x ab allows us to improve absorption but reduces the useful part of the computational domain. Values of τ ab are taken from the considerations of computational efficiency and minimization of the tails' influence on the bullet formation process. Thus, we can reduce the computational domain along the transverse coordinates. It is obvious that in such calculations the integrals (2.7)-(2.8) may not be conserved. Typically we take the domains L x = 20, L τ = 30 which are longer than those presented in Fig 1. It seems these dimensions are rather redundant. In fact, in the process of bullet trapping the pulse-beam spreads in time and x and parts of it which are close to the boundaries of absorption domains may be of intensities comparable with the peak ones. Chosen dimensions allow us to minimize this effect. At the same time they are shorter by an order of magnitude than those that should be taken in the absence of absorption layers.
Our computations are carried out on the basis of the difference scheme described in the Sec. 3. We choose also the following computational parameters: the step size along the propagation coordinate is h z = 0.001, while the step sizes along space and time coordinates are h x = 0.05, h τ = 0.05. Typical size of the computational domain along z axis in our simulation is L z = 500. Light bullets in a gradient waveguide

Simulation of light bullet trapping in different regimes
Since anomalous dispersion supports light bullets in 2D+1 media even without waveguide [6,8], our principal interest is to investigate waveguide at normal group velocity dispersion. The simulation and observation are conducted up to z = 500l nl . Such distance is an optimal one to investigate effects appearing due to the interference of diffraction, dispersion, nonlinearity and waveguide geometry.
In the first series of numerical experiments, we use a pulse of the Gaussian form at both frequencies (4.1) as a two-component input radiation. This case implies the phase-matching conditions k 2 = 2k 1 or Δk = 0. We have demonstrated that the regime of robust soliton propagation is established with the termination of energy exchange approximately at 200l nl . This propagation is accompanied by regular in-phase oscillations of peak intensities at both harmonics. The generalized phase F = 2φ 1 −φ 2 −Δkz oscillates near optimal value equal to −π. It also confirms the formation of a stable parametrically coupled pair of solitons and means that we observe the reactive regime of the classic parametric solitons when there is no energy exchange between the fundamental and second harmonics [5]. We have examined the temporal and spatial profiles at z = 500l nl as well. Due to waveguide influence the spatial profile is narrower than the temporal one. The calculated profiles have been compared with those of the Gaussian form having the same amplitudes, widths, and durations. This comparison demonstrates quite a good match between them. Thus, the soliton solution in this case has the form close to the Gaussian one. In total, in this computational series we have showed the trapping of a two-component "breathing" light bullet and its stability. One should underline that it forms at normal dispersion due to waveguide geometry only.
Consecutive processes of SHG and two-component light bullet formation are in the focus of the second series of our numerical simulation. In this experiment we launch a pulse in the Gaussian form (4.1) at the fundamental frequency only. Fig 2(A) demonstrates the dependence of the peak intensities of both harmonics on the propagation coordinate z. Firstly, energy transfer and second harmonic generation are observed. Approximately at a distance of 100l nl we see that both waves are trapped in a coupled pair and the regime of robust two-component soliton propagation is established. As in the previous case wave propagation is accompanied by regular in-phase oscillations of peak intensities at both harmonics (Fig 2(A)). Peak amplitude maxima are lower than in the previous case due to the lower initial total energy. Fig 2(B) shows the generalized phase oscillations near optimal value equal to π. In contrast to the previous case, the optimal generalized phase changes a sign. Note that the generalized phase is calculated as arccos function in the following way. If at a current point it is close to 2π, then at the next point the value 2πm is added to the calculated value in order to avoid a jump. Thus, at long distances the phase will be approximately 2πm (m>1 if it increases, or m<−1, if it decreases). Then, when a soliton is trapped, the phase will be close to π+2πm. The dimensionless temporal (Fig 2(C)) and spatial (Fig 2(D)) widths of the harmonics oscillates in-phase to each other and in-anti-phase to the intensity. Due to the influence of waveguide geometry spatial profile (Fig 2(F)) is narrower than a temporal one (Fig 2(E)). Fig 2(G) and Fig 2(H) show the evolution of temporal and spatial intensity profiles, correspondingly, at the fundamental frequency along the longitudinal coordinate z. Fig 3 illustrates the consecutive second harmonic generation and two-component light bullet trapping provided a small phase-mismatching. We see that the process of light bullet propagation is rather stable at Δk = 0.1 (Fig 3(A)). If we increase the value of the phase velocity detuning up to 0.5, it distorts the picture of the formation of a bullet due to a sufficient decay of the resulting pulse intensity. Even at Δk = 0.15 we observe a soliton close to the stability limit (Fig 3(B)). High values of phase mismatching prevent SHG. Fig 3(C) represents maximum amplitudes of both harmonics in time and transversal coordinate, which are averaged over the distance from z = 100 to z = 200. It can be seen that, starting from Δk = 0.2, the averaged maximum amplitudes drastically fall and soliton formation does not occur.
As it was discussed in the Sec. 2.1 the transverse susceptibility profiles may have different form. The next series of experiments deals with the Lorentzian (Fig 4), tangential ( Fig 5) and Gaussian (Fig 6) modes, respectively. The consecutive SHG and two-component light bullet formation is demonstrated in all three cases. The most evident confirmation of this fact can be received when we analyze the behavior of generalized phase-at the distance of 50l nl it begins to oscillate near optimal value equal to π. Since the waveguide in these cases effectively influences only on pulse-beams of a finite width, the soliton tails can gradually move away from the pulse center, a part of the pulse energy gradually goes to the periphery. This results in forming a quasi-stable bullet little by little losing its energy. In Fig 4(C) and Fig 5(C) it can be seen that the time profile also strongly deviates from the Gaussian one. At the same time, the generalized phase experiences significant oscillations, which indicate an obvious deviation of this regime from the soliton one. For the Gaussian profile (Fig 6), the deviations from the basic test case (Fig 2) are less noticeable. It may be explained by a good match between initial beam size and waveguide size. Fig 7 illustrates a remarkable property of the investigated system to capture both pulses in a waveguide and to trap a two-component light bullet while initially the waveguide presents at one frequency only. In this experiment we launch input pulses of the Gaussian form at both frequencies (4.1), provided waveguide with parabolic profile (2.3) is just at the fundamental frequency. We see that the bullet energy gradually decreases down to bullet disappearing.   However, one should note that the distance of total bullet destruction is not so small-it is equal approximately to 30 dispersion lengths. It is worthwhile to underline once again that provided normal dispersion, trapping of a two-component light bullet is impossible in the absence of a waveguide at least at the fundamental frequency.

Conclusions
To trap multi-component light bullets at normal dispersion of group velocities one should balance nonlinearity, dispersion, and diffraction. It is a challenging problem requiring the presence of an additional "power" compressing the pulses. Waveguides may play such a positive role, therefore, it was reasonable to investigate in detail a possibility of spatial-temporal solitons formation and propagation in a waveguide with quadratic nonlinearity.
Using the quasi-optical approach, we introduce a system of equations describing the propagation of pulse-beams in gradient waveguides. Numerical method for direct system simulation in the planar case is developed. A distinctive feature of the method is the preservation of the integrals of motion which are intrinsic to the governing system. Numerical simulation has been performed for various sets of parameters. Since light bullets at anomalous dispersion form and propagate even without waveguide, our principal interest was to investigate waveguide at normal group velocity dispersion. The formation of optical bullets has been shown when launching the input Gaussian pulse at both frequencies and at the fundamental frequency only as well. Besides that, we have computed the cases of phase mismatch and various types of waveguide profiles, in which soliton solutions also trapped. If a waveguide is only at the fundamental frequency, the formation of soliton-like solutions manifests itself too but its energy is gradually decreased. Nevertheless, as a whole it spreads to tens of dispersion lengths. In general, the simulations performed show a possibility of the existence of optical bullets and quasi-soliton solutions in the presence of normal dispersion in quadratically nonlinear waveguides under various conditions.