{'fraction_non_alphanumeric': 0.04636109167249825, 'fraction_numerical': 0.0233554933519944, 'mean_word_length': 4.488861148454005, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Universal fractal structures in the weak interaction of solitary waves in generalized nonlinear Schr\"{o}dinger equations

Weak interactions of solitary waves in the generalized nonlinear Schr\"{o}dinger equations are studied. It is first shown that these interactions exhibit similar fractal dependence on initial conditions for different nonlinearities. Then by using the Karpman-Solov'ev method, a universal system of dynamical equations is derived for the velocities, amplitudes, positions and phases of interacting solitary waves. These dynamical equations contain a single parameter, which accounts for the different forms of nonlinearity. When this parameter is zero, these dynamical equations are integrable, and the exact analytical solutions are derived. When this parameter is non-zero, the dynamical equations exhibit fractal structures which match those in the original wave equations both qualitatively and quantitatively. Thus the universal nature of fractal structures in the weak interaction of solitary waves is analytically established. The origin of these fractal structures is also explored. It is shown that these structures bifurcate from the initial conditions where the solutions of the integrable dynamical equations develop finite-time singularities. Based on this observation, an analytical criterion for the existence and locations of fractal structures is obtained. Lastly, these analytical results are applied to the generalized nonlinear Schr\"{o}dinger equations with various nonlinearities such as the saturable nonlinearity, and predictions on their weak interactions of solitary waves are made.


I. INTRODUCTION
Solitary wave interactions are a fascinating and important phenomenon for both physical and mathematical reasons. Physically, such interactions have arisen in a wide array of disciplines such as water waves [1], optics [2,3,4,5,6,7], and Josephson junctions [8]. For instance, in soliton-based fiber communication systems, optical pulses traveling in different frequency channels pass through each other, giving rise to collisions (strong interactions) of solitary waves. In the same frequency channel, neighboring optical pulses interfere with each other through overlapping tails, giving rise to weak interactions of solitary waves. Motivated by these physical applications, solitary wave interactions has been studied extensively in both the mathematical and physical communities. If the system is integrable, collisions of solitons are elastic [1], and weak interactions of solitons exhibit interesting yet simple behaviors [2,9,10,11,12]. However, in non-integrable systems, solitary wave interactions can be far more complex. The first sign of this complexity was reported by Ablowitz, et al. [13] for kink and anti-kink collisions in the φ 4 model where, inside the trapping interval, a reflection window was found. Later extensive numerical studies on this model by Campbell, et al. [14,15,16,17] revealed that in fact, sequences of two-and more-bounce reflection windows exist, and the physical mechanism for these refection windows is a resonant energy transfer between the translational motion and internal modes of kinks/antikinks. Anninos, et al. [18] pointed out further that there is a fractal structure in kink-antikink collisions. Using a collective-coordinate (i.e., variational) approach, they derived a set of fourthorder ordinary differential equations (ODEs) for these collisions, and these ODEs exhibit qualitatively similar fractal structures as in the φ 4 model (a comprehensive review on kink-antikink collisions in φ 4 -type equations can be found in [19]). These complex dynamics turn out to be not restricted to kink-antikink collisions. Indeed, similar phenomena have been reported on kink-defect collisions in the sine-Gordon and φ 4 models [20,21,22], as well as vector-soliton collisions in the coupled nonlinear Schrödinger (NLS) equations [23,24,25]. Furthermore, fractal scattering has also been reported on weak interactions of breathers in a weakly discrete sine-Gordon equation [26] and weak interactions of solitary waves in a weakly discrete NLS equation [27]. Recently, Goodman and Haberman [28,29,30] provided a deep analy-sis on the collective-coordinate models (ODEs) for kinkantikink collisions in the φ 4 model [18], kink-defect collisions in the sine-Gordon model [20], and vector-soliton collisions in the coupled NLS equations [25,31], using sophisticated dynamical system techniques. They derived analytical formulas for the locations of reflection-window sequences, which agree qualitatively with numerical results on the original partial differential equations (PDEs). Their results shed much light on the origins of these window sequences and fractal structures, especially from a mathematical point of view.
Despite the above progress on solitary wave interactions, our understanding on these phenomena is far from satisfactory. On the collision of solitary waves, the analysis done so far were all based on approximate collectivecoordinate approaches, hence the reduced ODE models can only provide qualitative results at best. Many features reported in the ODE models can not be seen in the PDE simulations, thus it is not possible to make a reliable prediction on the collision dynamics based on those ODE models and their analysis. In addition, the ODE models obtained from the collective-coordinate approaches not only are complicated, but also differ significantly from one PDE system to another. This forced previous researchers to analyze each PDE and its reduced ODE systems on an individual basis, which prevents an overall understanding on collision processes of solitary waves. On the weak interaction of solitary waves, the situation is even less satisfactory. The fractal nature of this weak interaction was reported only for systems which are weakly perturbed integrable systems (sine-Gordon and NLS equations, to be more specific) [26,27]. It is not known yet whether similar phenomena arise in strongly non-integrable equations. More seriously, the previous work on this subject is largely numerical. No analysis has been attempted yet (not even the approximate collectivecoordinate studies). Thus an analytical understanding on weak interactions of solitary waves is a completely open question.
In this paper, we study weak interactions of solitary waves in a whole class of generalized NLS equations (with arbitrary nonlinearities) both analytically and numerically. These generalized NLS equations are not weak perturbations of the NLS equation in general. First we show by direct PDE simulations that these weak interactions for different nonlinearities exhibit similar fractal structures on initial parameters of solitary waves. This establishes that fractal scattering is a common feature of weak interactions in this class of generalized NLS equations. Next, we rigorously derive a universal system of dynamical equations (ODEs) for the velocities, amplitudes, positions and phases of interacting solitary waves in this class of PDEs by the Karpman-Solov'ev method. This universal ODE system is remarkably simple, and it contains only a single parameter which depends on the individual PDEs (after variable rescalings). When this parameter is zero, these dynamical equations are integrable, and their exact analytical solutions are derived. When this parameter is non-zero, the dynamical equations are found to exhibit fractal structures for a wide range of initial conditions. These fractal structures match those in the original PDEs both qualitatively and quantitatively, thus the universal nature of fractal scattering in the weak interaction of solitary waves is analytically established. We further explore the origin of these fractal structures. Our numerical studies on the ODE system show that these fractal structures bifurcate from the initial conditions where the solutions of the integrable dynamical equations develop finite-time singularities. Based on this observation, we present an analytical criterion for the existence and locations of fractal structures. One corollary from this criterion is that when the initial separation velocity is above a certain threshold value, fractal structures should disappear -a prediction which agrees with our PDE numerics as well as previous numerics on the weakly discrete NLS equation (see Fig. 6 in Ref. [27]). Lastly, we apply these analytical results to the generalized NLS equations with various nonlinearities such as the cubic-quintic, exponential and saturable nonlinearities, and make detailed predictions on the dynamics of their weak interactions.
This paper is structured as follows. In Sec. 2, we describe individual solitary waves in the generalized NLS equations. In Sec. 3, we present direct PDE simulation results on weak interactions of solitary waves in the generalized NLS equations with two different nonlinearities, and reveal the common (universal) fractal structures in this class of PDEs. In Sec. 4, we analytically derive a universal system of dynamical equations (ODEs) for parameters of interacting solitary waves using asymptotic methods, and show that these ODEs accurately describe the weak interactions in the PDEs. In Sec. 5, we solve this ODE system analytically when the single parameter in this system is equal to zero (which is the integrable case). In addition, we derive explicit conditions for the solutions of the integrable ODE system to develop finitetime singularities. In Sec. 6, we show that fractal structures appear in this ODE system when its parameter is non-zero, and explore the origin of these fractal structures. In Sec. 7, we apply the analytical results to the generalized NLS equations with various nonlinearities. In Sec. 8, we summarize the results of the paper and make some further remarks.

II. PRELIMINARIES
The generalized NLS equation is where F (·) is a real-valued algebraic function with F (0) = 0. Equation (2.1) supports solitary waves of the form where Φ(θ) is a positive function which satisfies the following equation and β (> 0), V, x 0 , σ 0 are real constants. For convenience, we introduce the notations Physically, β is the propagation constant which is related to the solitary-wave amplitude (henceforth, we call β an amplitude parameter), φ is the phase of the solitary wave, σ 0 is its initial phase, ξ is its center position, V is its velocity, and x 0 is its initial position. The solitary wave is characterized uniquely by its four parameters: V, β, σ 0 and x 0 . The asymptotic behavior of this solution at infinity is where c is the tail coefficient which is determined by the nonlinear function F and propagation constant β. We define the power of the solitary wave as which plays an important role in the linear stability of the solitary wave. For general functions F , the analytical formulas for Φ, P and c are not available. But for some special nonlinearities, one can find the analytical solutions. For instance, for the cubic-quintic nonlinearity the analytical formulas for Φ, P and c are [6,32]  When two solitary waves are placed adjacent to each other, they would interfere through tail overlapping. In this case, the initial condition is where Φ satisfies Eq.(2.3). Here "0" in the subscript represents the initial value of the underlying parameter. For convenience, we assign the left solitary wave with index k = 1, and the right solitary wave with index k = 2.
To study the weak interaction between these two solitary waves, we require that the two solitary waves are both stable, well separated, and having almost the same velocities and amplitudes. Introducing notations the above requirements then amount to Here, P β ≡ dP/dβ > 0 corresponds to the Vakhitov-Kolokolov criterion for the linear stability of solitary waves in Eq. (2.1) [6,33]. Below, we numerically study the weak interaction of solitary waves in Eq. (2.1). This equation is numerically integrated by the pseudo-spectral method coupled with the fourth-order Runge-Kutta integration along the time direction. Since each solitary wave has four parameters, we have eight parameters in the initial conditions. Due to the phase, translation and Galilean invariances of Eq. (2.1), we can fix σ 0,1 = 0, x 0,1 + x 0,2 = 0 and V 0 ≡ (V 0,1 + V 0,2 )/2 = 0 without any loss of generality. Also, for simplicity, we take ∆V 0 = V 0,2 − V 0,1 = 0 in all our simulations of this section, i.e. the two solitary waves are initially at rest. This leaves four free parameters in the initial conditions (3.1): We also define the collision timet as the time when the two solitary waves are the closest (i.e. the separation distance between peaks of the two solitary waves the smallest) during interactions. The life time of interaction is defined as the time length from the beginning (t = 0) to the collision timet, which is equal to the collision time in value. Thust will be used to denote the life time as well. Of the four parameters in initial conditions, we will fix β 0 , ∆β 0 and ∆x 0 , and use ∆φ 0 as the control parameter and vary it continuously between 0 and 2π. At each  A. Weak interactions for the cubic-quintic nonlinearity Our first example of nonlinearity is the cubic-quintic nonlinearity (2.7), which arises in a wide array of physical systems such as optics [6] and boson condensates [32,34]. In this nonlinearity, we set α = 1, γ = 0.04. (3.5) It is easy to verify that all solitary waves (2.8) in this case are linearly stable using the Vakhitov-Kolokolov criterion. In our simulations, we set ∆x 0 = 10 and β 0 = 1. The x interval is 70 units wide, discretized by 512 grid points; and the time step size is 0.004. We first study the nonequal-amplitude case, and take ∆β 0 = −0.065. The ∆V ∞ versus ∆φ 0 diagram is shown in Fig.1. The prominent structures in this graph can be split into two regions: one region is −0.34 < ∆φ 0 < 2.5, and the other region is 2.9 < ∆φ 0 < 3.3. The structures in these two regions turn out to be quite similar (except a horizontal reflection with respect to a vertical axis), thus we focus on the larger region −0.34 < ∆φ 0 < 2.5 below. The main structure in this region forms a sequence of hills; their widths get smaller from the right to the left, and their heights are about the same. These hills will be called the primary hills. This primary-hill sequence converges to the accumulation point ∆φ 0c = −0.339. In order to see this hill sequence near the accumulation point ∆φ 0c more clearly, we zoom in the region [−0.35, 0.4], and the zoomed-in diagram is shown in Fig.2. In this figure, the cascading sequence can be seen very clearly (see Fig.2(a)). In Fig.2(b), the corresponding life-time diagram is displayed. We can see that on the same hill, interactions have roughly the same life time. On different hills, life times are different: hills closer to the accumulation point ∆φ 0c have longer life times. Between hills, even longer life times can be seen, suggesting more complex dynamics there. To explore differences in interaction dynamics on different hills, we select three points ∆φ 0 = 0.1759, −0.0057, −0.1053 (marked in Fig.2(a) by circles) on three adjacent primary hills. These points are at the same relative positions (roughly halfway between the peak and bottom) of the respective hills. At these points, the interaction dynamics is plotted in Fig.2(1-3).
Here only the separation distance ∆ξ versus time t graphs are shown. We find that these three dynamical processes are similar, except that the oscillation times before final separation differ by one from one hill to the next. The life timest n of interactions on this primary hill sequence are found to be an almost perfect linear function of the hill index n as ωt n = 2nπ + δ, (3.6) where the least-square linear fit gives Here the life time of each primary hill is measured numerically at the relative location of that hill shown in Fig.2(a) by circles. This life-time formula has the same form as those for all window sequences reported before [15,16,20,23,29]. In addition to the primary hill sequence as described above, Fig.1 also possesses higher-order structures between primary hills. To demonstrate, we first isolate the long interval [−0.35, 2.5] in Fig. 1 and re-plot that part of the graph in Fig.3(a). Then we zoom into its subinterval [0.91, 0.995], which is between the two largest primary hills in Fig.3(a). The zoomed-in graph is shown in Fig.3(b). We see that the zoomed-in graph is similar to Fig.3(a), but the cascading direction has reversed. This behavior is analogous to that reported in [18,24] for the φ 4 model and the coupled NLS equations. The main structure in this zoomed-in window is again two sequences of hills, accumulating to the left and right respectively. We call them secondary hills. Between secondary hills, we see even higher-order structures. To see these structures more clearly, we zoom into the sub-interval [0.9657, 0.96785], which is between the two largest secondary hills in Fig.3(b). The zoomed-in graph is shown in Fig.3(c). We see that it is again similar to Fig. 3(b) but with a reversed cascading direction. One can zoom into the regions between these tertiary hills in Fig.3(c) further, and will get even higher order structures which are similar to the ones shown in Fig.3. Thus Fig.3(a) is a fractal structure! We have also explored the interaction dynamics on this fractal. To demonstrate, we pick three ∆φ 0 values 1.511, 0.95248, 0.9669 which are at the same relative positions of the fractal (roughly halfway between the peak and bottom of the widest hills) in Fig.3(a)-(c) (marked by circles). The interaction dynamics at these three points are displayed in Fig.3(1-3) respectively. Here the positions of maximum amplitudes of the interacting waves are plotted against time. We see that these dynam- The exit velocity versus initial phase difference graph in the equal initial amplitude case of the cubic-quintic NLS equation 2.1), (2.7). The cubic and quintic nonlinearity coefficients as well as the initial conditions are the same as in Fig. 1, except that ∆β0 = 0 now. ical patterns are clearly similar, except that the numbers of oscillations before final separation are different.
In the above numerical simulations, the two solitary waves have different initial amplitudes (∆β 0 = −0.065). We have also studied interactions of equal-amplitude solitary waves, i.e. with ∆β 0 = 0, while keeping the other parameters the same. In this case, the graph of exit velocity ∆V ∞ versus initial phase difference ∆φ 0 is shown in Fig.4. This graph is symmetric with respect to ∆φ 0 for obvious reasons. Examination of this graph shows that it is also a fractal. Thus fractal dependence arises in weak interactions of both equal and non-equal amplitude solitary waves.

B. Weak interactions with exponential nonlinearity
To explore whether the above fractal structures for weak interactions persist or not with other types of nonlinearities, we consider in this subsection a different type of nonlinearity -the exponential nonlinearity, with Here, −1 is introduced into this function to meet the condition F (0) = 0. Note that this nonlinearity does not have any parameters. Throughout this subsection, we set the initial separation ∆x 0 = 8, and average propagation constant β 0 = 2.3. We study two cases, one for non-equal amplitudes with ∆β 0 = −0.045, and the other for equal amplitudes with ∆β 0 = 0. For both cases, the control parameter is ∆φ 0 as before. In our simulations, the x interval was 70 units wide, discretized by 512 grid points. The time step size was 0.002. The ∆V ∞ verse ∆φ 0 graphs for both cases are plotted in Fig. 5. We have verified that both graphs in this figure are fractals. Comparing these fractals with those in Figs. 1 and 4 of the cubic-quintic nonlinearity, we see that the fractal structures for these two different nonlinearities are very similar. The only major difference between them is in the non-equal amplitude case, where there is only one primary hill sequence (accumulating toward the left) for the exponential nonlinearity, while there are two primary hill sequences for the cubic-quintic nonlinearity. It is remarkable that two very different nonlinearities exhibit quite similar fractal dependence on initial conditions. Thus fractal scattering appears to be a universal feature in weak interactions of Eq. (2.1) rather than an accident. This leads us to the following questions: how can we analytically establish the universal nature of fractal scatterings for Eq. (2.1) with general nonlinearities? how can we analytically explain the major differences of fractals for different nonlinearities? These questions will be answered in the following sections.

IV. DYNAMICAL EQUATIONS
To study weak interactions analytically, we use the Karpman-Solov'ev method [9] by treating the interference as a small perturbation to each solitary wave (see also [35]). This method has been successfully used before on the NLS equation [2,9,10,35], the modified NLS equation [11], the Manakov equations [12], as well as the (non-integrable) coupled NLS equations [36]. To proceed, we first need to consider the evolution of a single solitary wave in the perturbed generalized NLS equation where function G is a perturbation term, and ǫ is a small parameter. Without perturbations (ǫ = 0), the solitary wave (2.2) is an exact solution of Eq.(4.1), and its internal parameters V, β, σ 0 , x 0 are time-independent. When the perturbation is turned on, these internal parameters of the solitary wave will evolve slowly on the time scale T = ǫt. The multiple-scale perturbation theory for this slow evolution is well known [36,37]. We write the perturbed solution as Here V (T ), β(T ), σ 0 (T ), x 0 (T ) are all functions of slow time T . Next, we will derive the dynamical equations (ODEs) for the slow-time evolution of these parameters. Substituting (4.2) into (4.1), we get the equation forΦ as where φ is defined in Eq. (2.4). We expand the amplitude functionΦ into a perturbation serieŝ The equation at order ǫ 0 is satisfied automatically since Φ satisfies Eq.(2.3). At order ǫ, the equation forΦ can be written as (4.8) and Here the superscript "*" represents complex conjugation. Operator L has two eigenfunctions and two generalized eigenfunctions associated with the zero eigenvalue, with the relations In order for the inhomogeneous solution Ψ of the firstorder equation (4.6) to be non-secular at large time, the inhomogeneous term in Eq.(4.6) must be orthogonal to the above eigenfunctions and generalized eigenfunctions of the zero eigenvalue, i.e., under the inner product defined as Here F † k is the Hermitian of F k . Evaluating the four integrals in Eq.(4.12), the slow-time evolution equations for parameters V (T ), β(T ), σ 0 (T ), x 0 (T ) will be obtained. These evolution equations can be written as These equations will be critical for the development of weak interaction theory of solitary waves below. Now, we consider the weak interaction of two solitary waves. Here the tail overlapping can be considered as a small perturbation which causes the internal parameters of each solitary wave to evolve on a slow time scale ǫt. Here ǫ is the magnitude of tail overlapping which is exponentially small with solitary wave spacing ∆ξ. We will not introduce ǫ explicitly in the next analysis. To the leading order, the interacting solution is simply a superposition of two solitary waves, where all parameters slowly vary over time. Picking up the dominant interference terms, each solitary wave is governed by the following perturbed generalized NLS equations: where In this paper, we only study the weak interaction, so conditions (3.4) are assumed. Since |∆V | ≪ 1, the phase difference which is independent of θ. Now We apply the above solitary wave perturbation theory to Eq. (4.19). In this problem, Substituting (4.22) into Eqs.(4.14)-(4.17), we obtain the following dynamical equations where P k , k = 1, 2 are powers of the two individual solitary waves. These equations can be simplified greatly. Due to assumptions (3.4), and noticing that Φ(θ) and Φ β (θ) are even functions of θ, the leading-order terms of the above integrals can be explicitly calculated. For instance, With the above simplifications, the dynamical equations reduce to where P is the power of the solitary wave with propagation constant β, and D 1 , D 2 are defined in (4.29),(4.30).
From the above equations, we find that and (4.43) Physically, this latter symmetry corresponds to the interchange of the left and right solitary waves in the PDE evolutions, which of course does not change the interaction outcome. Eq. (4.42) also has the property that if ψ(τ ) is a solution, so is ψ(τ ) + 2nπ for any integer of n. This reflects the fact that in the PDE system, solution evolution remains the same if the phase difference between the solitary waves changes by a multiple of 2π. The dynamical equations (4.42) are asymptotically accurate in describing weak interactions in the PDE system (to the leading order) when the spacing ∆ξ is large. Surprisingly, even when the two solitary waves come close to each other, Eq. (4.42) can still describe the interaction process very well. This is analogous to weak interactions in the (integrable) NLS equation [2,9]. Below, we make detailed comparisons between the ODE solutions of Eq. respectively. In both cases, the initial phase difference ψ 0 is the control parameter as in PDE simulations.  Next we will turn our attention to the ODE system (4.42), and analyze its solution dynamics in more detail. In particular, we would like to understand why fractal structures arise in this system, and how to analytically predict their locations and other main features. 2) The general solution of this equation is where and The asymptotic behaviors of these solutions as τ → ∞ can be easily determined. Let where a, b, c, d are real constants, then the following leading-order asymptotic expressions for the solution can be obtained when τ → ∞: From these asymptotic expressions, we see that when a = 0, the two solitary waves eventually move away from each other with exit velocity 2|a|; when a = 0 but b = 0, the solution is time-periodic for both d = 0 and d = 0, the difference being that in the former case, the periodic solution exhibits finite-time singularities (where ζ = Re(Y ) → ∞), while in the latter case, the solution has no singularities; when a = b = 0, the two solitary waves eventually separate logarithmically, and the exit velocity is zero. As an example, we take the initial conditions (4.44). In this case, the graph of exit velocity −ζ ∞ (= 2|a|) versus ψ 0 is plotted in Fig. 9 (bottom panel). This graph is smooth everywhere, except at ψ 0 = 0, ±π where it has a cusp (due to the absolutevalue function in |a|). The squares and diamonds on this graph will be explained later. Clearly, this graph has no fractal structure anywhere. Thus, fractal dependence is a signature of the dynamical system (4.42) when it is nonintegrable (with ε = 0), not when it is integrable (with ε = 0).
The above asymptotic states do not tell the full story about the solution dynamics in the integrable system. For instance, for the case of a = 0, even though the solution has a benign-looking asymptotics (5.8) as τ → ∞, the solution can still develop a singularity (where the separation ζ → ∞) at a finite time. These solutions with finite-time singularities turn out to be critical for the appearance of fractal structures in the non-integrable system, as our numerics in the next section will indicate. Thus we analyze these singularity solutions in more details below. The necessary and sufficient conditions for singularities in solution (5.3) are that cosh(C 1τ + C 2 ) = 0, (5.12) andτ > 0, whereτ is the time of singularity. Ifτ < 0, i.e., singularities in the solution occur at a negative time, such singularities are irrelevant for the time evolution of Eq. (4.42) and need not be considered. The solutions of Eq. (5.12) are C 1τ + C 2 = 1 2 (2n + 1)πi, n = 0, ±1, ±2, · · · . (5.13) This is a complex-valued relation, which gives two real relations onτ , C 1 and C 2 . When a = 0, i.e., C 1 is not purely imaginary, we find by separating the real and imaginary parts of Eq. (5.13) that the solution (5.3) has a single finite-time singularity of the type ln(τ −τ ) if the following conditions are satisfied: S ≡ Im(C * 1 C 2 ) Re(C 1 ) = 1 2 (2n + 1) π, n = 0, ±1, ±2, · · · , (5.14) Here Re(·) and Im(·) represent the real and imaginary parts of a complex number. When a = 0 (b = 0), singularity solutions exist if d = 0. These solutions have an infinite number of finite-time singularities of the type ln(τ −τ ), as the formula (5.9) indicates. Physically, at the time of singularityτ , the two solitary waves strongly collide, thusτ is the collision time. Whether conditions (5.14) and (5.15) can be satisfied depends on the initial conditions (which determine the C 1 and C 2 values, see (5.4), (5.5)). In the text below, we will call initial conditions (ζ 0 ,ζ 0 , ψ 0 ,ψ 0 ) which satisfy Eqs. (5.14) and (5.15) as singularity points. At singularity points, solutions of the integrable dynamical system (4.42) develop finite-time singularities.
To demonstrate how to determine singularity points in the initial-condition space, we take initial conditions (4.44) of Fig. 6(a) as an example. Here ψ 0 is a control parameter. With these initial conditions, the graph of function S(ψ 0 ) is plotted in Fig. 9 (top panel). This graph has a maximum 0.96. As ψ 0 → 0 + or π − , S(ψ 0 ) → −∞. As we can see from this graph, for any value of n ≤ −1, Eq. (5.14) has two roots, ψ (1) 0,n and ψ (2) 0,n . We have checked that these roots satisfy the other singularity condition (5.15), thus these ψ (1) 0,n and ψ (2) 0,n values are singularity points. It is noted that the graph of function S(ψ 0 ) also has another piece in the interval π < ψ 0 < 2π, which is the mirror image of that shown in Fig. 9 around the point ψ 0 = π. But in that interval,τ < 0, not satisfying the second singularity condition (5.15), thus we did not plot that piece of the graph in Fig. 9.
The above detailed analysis on singularity points was performed for the particular initial conditions (4.44) where the two solitary waves are initially stationary (ζ 0 = 0). What will happen ifζ 0 = 0? To answer this question, we fix ζ 0 andψ 0 as in Eq. (4.44), varyζ 0 , and examine how singularity points move in the (ψ 0 ,ζ 0 ) plane. The results are shown in Fig.10. The top curve corresponds to n = −1, the next curve corresponding to n = −2, and so on. All curves are bounded from both above and below except the top one (with n = −1). When n → −∞, these curves approach the accumulation curve plotted by the dashed line in Fig.10. Below this accumulation curve, there are no singularity points. The analytical formula for this accumulation curve can be easily derived. On this accumulation curve, C 1 must be pure imaginary, thus M =ζ 0cψ0 − e ζ0 sin ψ 0 = 0, (5.19) For the ζ 0 andψ 0 values of Fig. 10, we getζ 0c,min = −0.00389. Ifζ 0 <ζ 0c,min , there are no singularity solutions for any value of ψ 0 . Whenζ 0c,min <ζ 0 <ζ 0c,max , two infinite sequences of singularity points can be found. Whenζ 0 >ζ 0c,max , however, the number of singularity points becomes finite; this number gradually decreases (down to one) asζ 0 increases. The above calculations of singularity points and their accumulation curves were made for special choices of initial conditions ζ 0 = −10 andψ 0 = −0.01167 (see Fig.  10). In view of the importance of singularity points for fractal structures which we will reveal in the next section, we would like to discuss singularity points and their accumulation curves further for general initial conditions below.
First, we examine the accumulation curve in the (ψ 0 , ζ 0 ) plane for general initial conditions ζ 0 andψ 0 . In this general case, the accumulation curve (if it exists) is necessarily given by Eq.
Another general property about singularity points is that, at sufficiently large values ofζ 0 , there is a unique singularity point in the ψ 0 interval. The proof is as follows. It is easy to check that whenζ 2 0 +ψ 2 0 ≫ e ζ0 and |ζ 0 | > ∼ |ψ 0 |, functions S andτ have the following leadingorder asymptotic expressions, (5.30) and the relative errors are O e ζ0 /(ζ 2 0 +ψ 2 0 ) . From Eq. (5.29), we see that the rise of S value over the interval 0 ≤ ψ 0 ≤ 2π is π, which guarantees that Eq. (5.14) has a single solution in the ψ 0 interval for a single value of n. From Eq. (5.28), we see that whenζ 0 > 0 is sufficiently large,τ > 0 over the entire ψ 0 interval. Thus singularity conditions (5.14) and (5.15) admit a unique singularity point. We note by passing that whenζ 0 is sufficiently large negative,τ < 0 over the entire ψ 0 interval, thus there can not be any singularity points. This is consistent with the previous general property proved above.
It is noted that in the above classifications, case (2) is just a special case of case (3), and can be readily deduced from (3). Case (1) can be deduced from (3) as well under the limitψ 0 → 0. But cases (1) and (2) are important special cases, hence we listed them out separately.

VI. ORIGINS OF FRACTAL STRUCTURES IN THE NON-INTEGRABLE DYNAMICAL EQUATIONS
We have known from Figs. 6, 7 and 8 that the nonintegrable ODE system (4.42) exhibits hill sequences and fractal structures which coincide with those in the PDE simulations, but such structures do not exist when this ODE system becomes integrable. The natural question then is: where do the fractal structures in the nonintegrable system (4.42) come from? In this section, we will establish through careful numerics that these fractal structures bifurcate from singularity points of the integrable system.
To determine the origin of these fractals, we take the same initial conditions (4.44) as in Fig. 6(a), but gradually decrease the value of ε from 0.13665 of Fig. 6(a)  down to zero (the integrable case), then down further to negative values. In this process, we closely monitor how the fractal structure of Fig. 6(a) changes as ε decreases. The result is shown in Fig. 11. Here, the −ζ ∞ verse ψ 0 graphs are plotted at six decreasing ε values: ε = 0.13665, 0.036, 0.0036, 0, −0.0036 and −0.036. We see that as ε decreases from 0.13665 but above zero, primary hill sequences and the fractal regions between them persist and are clearly visible in Fig. 11(1, 2, 3). Indeed, we have zoomed into the sensitive regions between primary hills in each of Figs. 11(1, 2, 3), and obtained higher order structures which look very similar to those shown in Fig. 8. As ε → 0 + , our key observation is that, the peaks of individual primary hills as well as the nearby fractal regions collapse to sequences of points on the smooth −ζ ∞ curve of the integrable system (see Figs. 11(3, 4)). Closer examination tells us that, these sequences of points in Fig. 11(4) are nothing but the two sequences of singularity points of the integrable system which we plotted in Fig. 9! In other words, hill sequences and fractal structures in the non-integrable system bifurcate from the singularity points of the integrable system! However, this bifurcation is one-sided: as ε decreases below zero, no fractal regions appear, see Fig. 11 (5). A finite number of primary hills, reminiscent of primary hill sequences for positive ε values, do exist. But the whole graph is smooth, and it has no fractal structures inside (even the spike-looking parts of the graph in Fig. 11 (5) turn out to be smooth upon closer examination). Furthermore, as ε decreases further below zero, the number of primary hills keeps decreasing, and the graph becomes more smooth, see Fig. 11 (6). Thus, fractal structures are a signature of the non-integrable system (4.42) only for positive values of ε, not negative values of ε.
To further substantiate our claim on fractal structures of the non-integrable system bifurcating from singularity points of the integrable system, we tune initial conditions so that singularity points in the integrable system gradually disappear in the ψ 0 interval, and check if fractal structures in the non-integrable system disappear as well (for small ε). Specifically, we fix the ζ 0 anḋ ψ 0 values in Eq. (4.44) and tune theζ 0 value, as we did in Fig. 10. The ε value in Eq. (4.42) is taken as ε = 0.0036, which is very small. Thus, the non-integrable system is weakly perturbed from the integrable one. For the above initial conditions, singularity points of the integrable system have been displayed in Fig. 10 in the (ψ 0 ,ζ 0 ) plane. We gradually decrease theζ 0 value. For eachζ 0 , we numerically compute the exit velocity versus ψ 0 graph in the perturbed (non-integrable) system (4.42), and compare how this graph relates to singularity points of the integrable system in Fig. 10. To illustrate, we pick seven representativeζ 0 values, which are 0.00707, 0.00548, 0.00495, 0.00350, 0, −0.00350 and −0.00424 in decreasing order. These sevenζ 0 values are marked by horizontal dashed lines in Fig. 10. As we can see from that figure, at these sevenζ 0 values, the numbers of singularity points in the ψ 0 interval are 1, 3, 5, ∞, ∞, ∞, and 0 respectively. For each of these sevenζ 0 values, the corresponding exit velocity versus ψ 0 graph in the perturbed system (4.42) is shown in Fig. 12. We notice from this figure that the numbers of primary hills and fractal regions near these hills at theseζ 0 values are equal to 1, 3, 5, ∞, ∞, ∞, and 0 respectively -exactly like singularity points in the integrable system! In particular, when singularity points in the integrable system disappear, so do primary hills and fractal structures in the weakly perturbed non-integrable system. Furthermore, the locations of primary hills and fractal regions closely follow those of singularity points of the integrable system. Thus, the connections between them are unmistakable. Fig. 12, together with Fig. 11, establishes beyond doubt that primary hills and fractal structures in the non-integrable system (4.42) indeed bifurcate from singularity points of the integrable system.
The bifurcation of fractal structures from singularity points of the integrable system indicates that near such points, the solutions of the perturbed system (4.42) are very sensitive to initial conditions. To shed light on why this sensitivity occurs, we present some numerical results below. First, we look at the integrable system (with ε = 0). Taking the initial conditions as (4.44), evolutions of ψ versus τ at the singularity point ψ 0 = 0.98325 (marked in Fig. 9, bottom panel) and its left and right near neighbors ψ 0 = 0.92 and 1.05 are plotted in Fig.  13(a) using the solution formula (5.3). An interesting feature about these evolutions is that for initial ψ 0 values at the two sides of the singularity point, the phase functions ψ(τ ) have drastically different trajectories as they go through the time τ ≈ 700 where the two solitary waves interact strongly (this time is the singularity time of the singular solution at ψ 0 = 0.98325). For ψ 0 below the singularity point, the phase sharply (but continuously) decreases by 2π, while for ψ 0 above the singularity point, the phase sharply (but continuously) increases by 2π. Recall thatψ ∝ ∆β and it determines the relative energy (amplitude) distributions among the two solitary waves [see Eqs. (4.37), (4.40) and (4.41)], we know that on the two sides of the singularity point, the energies have opposite distributions among the two solitary waves during their strong interactions. However, after the interaction is completed, the asymptotic slopes of the three ψ(τ ) trajectories in Fig. 13(a) are almost the same, signaling that the interaction outcome is actually insensitive to the ψ 0 values (the roughly 2π difference between these phase trajectories does not affect the physical solutions). This is why outcomes of weak interactions in the integrable system (4.42) do not exhibit sensitive dependence on initial conditions (see Fig. 9 bottom panel). However, when system (4.42) is positively perturbed, the results are completely different. To demonstrate, we take ε = 0.0036 now, while the initial conditions (4.44) remain the same. In this slightly perturbed system, the solution develops finite-time singularity at the singularity point ψ 0 = 0.9695, which is the counterpart of the singularity point mentioned above in the integrable system. The phase function at this singularity point is plotted in Fig.  13(b) (solid line). (It is noted that in this perturbed case, we do not have exact solution formulas, hence this solution was obtained by numerically integrating Eq. (4.42). Due to the finite-time singularity in the solution, our numerical integration can not go beyond the singularity time τ ≈ 700. The solution beyond the singularity time, shown in Fig. 13(b) as dotted lines, was inferred from our numerics at nearby ψ 0 values.) On the two sides of the singularity point, we select two nearby values ψ 0 = 0.92 and 1.01. The phase trajectories at these φ 0 values are also plotted in Fig. 13(b). We see that as these trajectories go through the time τ ≈ 700, one sharply decreases by 2π, while the other sharply increases by 2π, similar to what happens in the integrable case (see Fig. 13(a)). However, after these sharp decreases/increases, the trajectories turn around and start to move in the opposite direction. Eventually, these trajectories approach drastically different asymptotic slopes (one positive and the other one negative in fact), indicating that the interaction outcomes are very different for these slight changes in the ψ 0 values. This is the phenomenon of sensitive dependence on initial conditions which occurs in the perturbed system (4.42) (with ε > 0), but not the integrable system (with ε = 0).
It is also enlightening to look at this sensitive depen-dence on initial conditions from the viewpoint of PDE evolutions. To illustrate, we take the cubic-quintic nonlinearity (2.7) in Eq. (2.1), and take α = 1, β 0 = 1. Then for the ε values and initial conditions used in the ODE simulations of Fig. 13, and in view of the variable rescalings (4.40) and (4.41), the corresponding PDE parameters for Fig. 13(a) (the integrable case) are γ = 0, ∆β 0 = −0.066016, ∆V 0 = 0, ∆x 0 = 10, and the corresponding PDE parameters for Fig. 13(b) (the perturbed case) are γ = 0.0010, ∆β 0 = −0.066, ∆V 0 = 0, ∆x 0 = 10. For these PDE parameters, the PDE evolution results (in the form of contour plots) at three ∆φ 0 values corresponding to those in the ODE simulations of Fig. 13 are displayed in Fig. 14. In the integrable (NLS) case (top row of Fig. 14), we take the three ∆φ 0 values exactly the same as those in Fig. 13(a), i.e. ∆φ 0 = 0.92, 0.98325, 1.05. In this case, at the lower ∆φ 0 value, the left solitary wave retains its higher energy at the collision time; at the singularity point of ∆φ 0 , the two waves completely coalesce at the collision time, signaling the singularity formation in the ODE system; at the higher ∆φ 0 value, the right solitary wave gets higher energy at the collision time. However, after interaction, the two waves always separate, and the right wave always gets higher energy, in all three cases. Recall that before interaction, the left wave has higher energy, thus we can call these interaction outcomes transmission. In these interactions, even though the intermediate process (especially the collision segment) rather strongly depends on the initial phase difference ∆φ 0 , the interaction outcome is insensitive to it. These PDE evolution results completely resemble the ODE simulations in Fig. 13(a). In the perturbed (non-integrable) case, the PDE simulation results are quite different from the integrable ones (as in the ODE simulations). In the perturbed case, we take the three ∆φ 0 values to be 0.92, 0.972 and 1.01. Notice that the first and third of these ∆φ 0 values are exactly the same as those in the ODE simulations of Fig. 13(b), while the middle ∆φ 0 value of 0.972 is slightly different from the corresponding ODE value of 0.9695. This slight difference in the middle ∆φ 0 value is necessary in order for the corresponding ODE and PDE simulations to exhibit the same behaviors, and this difference is due to the modeling error of the PDE evolutions by the ODE system (4.42). At ∆φ 0 = 0.92 (see Fig. 13(1)), the interaction outcome is similar to the integrable ones (top row of Fig.  13) in that it is also transmission. But at ∆φ 0 = 0.972 ( Fig. 13(2)), the two waves strongly coalesce, then form an oscillating bound state. At ∆φ 0 = 1.01, on the other hand, the two exiting waves have opposite energy distributions from Fig. 13(1); this interaction outcome can be called reflection. Thus, in the perturbed case, the interaction outcome is sensitive to initial conditions, which distinctively contrasts the integrable case. Again, these PDE evolution results for the perturbed case completely resemble the ODE simulations in Fig. 13(b).
The above ODE and PDE simulations corroborate the fact that the source of this sensitive dependence on initial conditions in the perturbed system lies in the finite-time singularities of solutions in the integrable dynamical system (4.42). From the PDE point of view, the origin of this sensitive dependence can be traced to the coalescing of the two solitary waves in the integrable PDE system. At the moment, our understanding on this sensitive dependence and fractal structures in the perturbed system is still very limited. For instance, we can not yet explain any quantitative details inside these fractal structures, nor can we explain why this sensitive dependence occurs only for one-sided perturbations of the integrable system (with ε > 0). These are non-trivial questions which merit further analysis, but they are beyond the scope of the present article.
The fact of primary hills and fractal structures in the non-integrable system (4.42) bifurcating from singularity points of the integrable system has far reaching consequences. One important consequence is that, the main features of primary hill sequences shown in Figs. 2(a) and 7(a) for PDEs and ODEs can now be analytically explained. For instance, the life-time formula (3.6) for primary hill sequences in the weakly perturbed system (4.42) is nothing but the analogous collision-time (singularity time) formula (5.16) for sequences of singularity points in the integrable system. To make a quantitative comparison between these formulas, we take the initial conditions (4.44) which was used in the PDE and ODE simulations of Figs. 2(a) and 7(a). When the time rescaling (4.41) is recovered, the collision-time formula (5.16) of the integrable system becomes 0.0839t n = 2nπ + π, (6.1) which compares very favorably with the life-time formulae (3.6), (3.7) and (4.46) in direct PDE and ODE simulations. The small differences in the ω and δ values between the analytical formula (5.16) and the PDE/ODE ones (3.6) are caused by the not-so-small value of ε = 0.13665. As ε → 0, these quantitative differences will vanish. Regarding the locations of individual hills in the primaryhill sequence, they are described by the formula (5.17) for singularity-point locations of the integrable system when ε ≪ 1. Note that the form of this formula is different from all previous ones on window sequences in solitarywave collisions [15,16,20,25,29]. For ε > 0, each primary hill is paired with a sensitive (fractal) region at its foot (see Figs. 11 and 12). Similar to primary hill sequences, the locations of these fractal regions are described by the same formula (5.17) in the limit ε → 0 + . The fact of primary hills and fractal structures bifurcating from singularity points of the integrable system also explains major features of interaction results in Fig.  5(a) for the exponential nonlinearity (3.8). We have noticed that, unlike Fig. 1, this graph has only one infinite sequence of primary hills accumulating toward the left (the right sequence of Fig. 1 is absent). This phenomenon is due to the fact that for the choices of initial conditions for Fig. 5(a), there is only one infinite sequence of singularity points in the integrable system. To see it, we first calculate the f, g and ε values for Notice thatψ 2 0 < 2e ζ0 , thus according to the classifications of singularity points in the end of the previous section (case 2), the integrable equation (4.42) with the above initial conditions has only one sequence of singularity points in the ψ 0 interval, accumulating to the left toward ψ 0 = 0 + . This is in perfect agreement with the primary-hill sequence of Fig. 5(a) from direct PDE simulations.
In many of the interaction results presented in this paper, the exit velocity versus ψ 0 graphs have infinite sequences of primary hills (see Figs. 1 and 5 for instance); when zooming into the sensitive regions between primary hills, one gets infinite sequences of secondary hills. It is important to understand that these two infinite sequences are pure coincidence, and are totally un-related. Each primary hill corresponds to a particular singularity point of the integrable system, thus the number of primary hills is equal to the number of singularity points in the integrable system. This number can be either infinite or finite, depending on the choices of initial conditions. For instance, Figs. 12(1, 2, 3) have 1, 3 and 5 primary hills, corresponding to the same numbers of singularity points on the top three dashed lines of Fig. 10. On the other hand, at the foot of each primary hill, there is always an infinite sequence of secondary hills (when ε > 0). In other words, secondary hills always exist as an infinite, not finite, sequence. For example, if one zooms into each of the three sensitive regions at the foot of the three primary hills in Fig. 12(2), one always gets an infinite sequence of secondary hills. Thus secondary-hill structures are un-related to primary-hill structures. If we zoom into the sensitive regions between secondary hills, we always get infinite sequences of tertiary hills which are very similar to the sequences of secondary hills both qualitatively and quantitatively (see Figs. 3(b, c)). This process can continue indefinitely. Thus, our conclusion is that sensitive regions between primary hills are fractal structures (in the sense that portions of these structures, when amplified, are the same as the strcutures themselves). But the whole graph with primary hills is not a fractal.

VII. APPLICATIONS TO THE GENERALIZED NLS EQUATIONS WITH VARIOUS NONLINEARITIES
In previous sections, we have shown that for the cubicquintic and exponential nonlinearities at selected parameters (α = 1, γ = 0.04, β 0 = 1 for the former, and β 0 = 2.3 for the latter), weak interactions of solitary waves exhibit hill sequences and fractal structures for a wide range of initial conditions, and the reduced ODE model (4.42) accurately captures these interaction dynamics both qualitatively and quantitatively. In this section, we consider a larger question: for a given form of nonlinearity in the PDE (2.1), can it exhibit fractal structures? For example, with the cubic-quintic nonlinearity (2.7), for what parameters α and γ can one possibly find fractal structures? This question can be answered by applying our previous results on the ODE model (4.42). For demonstration purpose, we will do so for three forms of nonlinearity: cubic-quintic, exponential, and saturable nonlinearities.
From the analysis of the ODE system (4.42) in the previous section, we have found that fractal structures in weak interactions can only occur for ε > 0, not for ε < 0. Thus, once we have obtained the functional dependence of ε on system parameters, it will quickly become clear when fractal structures can arise. The analytical expression for ε is given in Eq. (4.41). Notice that due to the Vakhitov-Kolokolov stability criterion [6,33], the solitary wave is linearly stable only when P β > 0, i.e. ε > −1. Below, we will use the ε formula in (4.41) to calculate ε for general system parameters in the three nonlinearities mentioned above.
First we consider the cubic-quintic nonlinearity (2.7). When α < 0, we found that P β is always negative, i.e. the solitary wave is always linearly unstable. Thus we only consider the α > 0 case below. In this case, it is easy to see from Eqs. (2.3) and 2.7) that by a rescaling of variables, we can make α = β 0 = 1. Thus the only remaining parameter for this nonlinearity is γ. Using the analytical formula (2.9) for P , we can obtain the dependence of ε on γ, which is plotted in Fig. 15(1). From this graph, we see that ε > 0 when γ > 0, and ε < 0 when γ < 0. Thus fractal structures in this cubicquintic model can appear only when γ > 0, not when γ < 0. If γ = 0, this cubic-quintic model reduces to the integrable NLS equation, and the dynamical equations (4.42) reduce to the integrable case (with ε = 0) studied in Sec. 5. In this integrable case, there is of course no fractal dependence in solitary wave interactions.
Next, we consider the saturable nonlinearity, which is common in optics (for instance, in photorefractive crystals [38]). Here one is added in the above formula to make F (0) = 0 (this does not affect the solitary waves and their interaction dynamics). In this case, ε also depends only on the propagation constant β. This dependence is computed numerically and plotted in Fig.15(c). We find that ε is negative for all values of β, thus fractal structures can not exist in weak interactions of solitary waves for this saturable nonlinearity. This conclusion is consistent with our earlier results for the cubic-quintic nonlinearity, as the saturable nonlinearity (7.1) resembles the cubic-quintic nonlinearity (2.7) with α > 0 and γ < 0. It is noted, however, that for the saturable nonlinearity, weak interactions of solitary waves can still exhibit some interesting structures as shown in Figs. 11 (5,6), but these structures are not fractal structures. From the above three examples (as well as the previous section), we see that the reduced ODE model (4.42) enables us to accurately predict when and where fractal structures and hill sequences appear in the space of initial parameters of solitary waves. Based on this reduced model, a global and universal understanding on weak interactions of solitary waves has been achieved for the generalized NLS equations (2.1) with arbitrary forms of nonlinearity.

VIII. CONCLUSION AND DISCUSSION
In this paper, we have analyzed weak interactions of solitary waves in the generalized nonlinear Schrödinger equations with general forms of nonlinearity. We have shown that these interactions exhibit similar fractal dependence on initial conditions for different nonlinearities.
To analytically explain these universal fractal structures, we derived a set of fourth-order dynamical equations for the solitary-wave parameters using asymptotic methods. A remarkable feature of these dynamical equations is that they contain only one parameter, which is dependent on the specific form of nonlinearity. When this parameter is zero, these dynamical equations are integrable, and the exact analytical solutions are derived. When this parameter is non-zero, the dynamical equations exhibit fractal structures which match those in the original PDEs both qualitatively and quantitatively. We have also investigated the origin of these fractal structures, and found that they bifurcate from the singularity points (i.e. initial conditions for singularity solutions) in the integrable system. Based on this observation, an analytical criterion for the existence and locations of fractal structures is obtained. Lastly, we applied these analytical results to the generalized nonlinear Schrödinger equations with various nonlinearities such as the cubic-quintic, exponential and saturable nonlinearities, and predictions on their weak interactions of solitary waves are presented.
Regarding the bifurcation of fractal structures from the integrable dynamical equations, even though we have established that this bifurcation occurs at the singularity points of the integrable system, more challenging questions are to comprehensively analyze how this bifurcation takes place, and to quantitatively predict the detailed geometric structures inside these fractals. This has not been done yet. Recently, Goodman and Haberman analyzed the approximate ODE models for the collisions of solitary waves in three physical systems where window sequences and fractal structures have been reported [28,29,30]. They found that the origin of window sequences and fractal structures in these systems lies in the crossing of the separatrix (homoclinic orbit). Analytical predictions on the locations of window sequences in the ODE models were derived as well. It is not clear at the moment whether similar analysis can be performed for our system (4.42). This question is beyond the scope of the present article, and will be left for future studies.