Chaoticons described by nonlocal nonlinear Schrödinger equation

It is shown that the unstable evolutions of the Hermite-Gauss-type stationary solutions for the nonlocal nonlinear Schrödinger equation with the exponential-decay response function can evolve into chaotic states. This new kind of entities are referred to as chaoticons because they exhibit not only chaotic properties (with positive Lyapunov exponents and spatial decoherence) but also soliton-like properties (with invariant statistic width and interaction of quasi-elastic collisions).

where the real positive function R(x) is the (nonlinear) response function, which must be symmetry for the existence of the soliton-like solutions 30 . The Hamiltonian of Eq. (1) 2 2 is conserved 29 . The nonlocal nonlinearity (the convolution integral) in Eq. (1) means that the wave-induced "potential" at a certain spatial point x, 2 , is determined not only by the wave q(x, t) at that point but also by the wave in its vicinity. This kind of nonlocal nonlinear response has been found in several systems, such as Bose-Einstein condensates 31,32 , atomic vapors 33 , nematic liquid crystals 15,16 , thermal susceptibilities 34 , etc. The stronger the nonlocality, the more extended the wave distribution contributing to the "potential" V 13,14,35 . Different from the (local) nonlinear Schrödinger equation [when R(x) = δ(x) in Eq. (1)], nonlocality has profound effects on the dynamics of solitons. For example, the interaction of two nonlocal solitons can have both a long-range mode 22,23,35 and a short-range mode 24,25,35 , but two local solitons interact with each other only in a short-range one 7,35 ; and the NNLSE [Eq. (1)] can support the multi-hump solitons with the Hermite-Gauss-type (HGT) profiles [18][19][20][21] , but the nonlinear Schrödinger equation admits only the single-hump solitons 5 . However, the NNLSE may not guarantee the existence of all high order HGT-solitons. The response function also plays an important role. The NNLSE with the Gaussian response function can support the HGT-solitons without upper threshold of the hump-number [18][19][20] . Contrastively, the NNLSE with the exponential-decay response function 13 only admits of the HGT-solitons with the hump-number less than five 20 . The crucial difference between such two kinds of response functions is 11,35 that the former is non-singular and the potential V can be simplified to a quadratic form in the limit of strong nonlocality, while the latter that can describe physically real materials is singular and the corresponding NNLSE cannot be generally reduced to a linear Snyder-Mitchell mode 12 .
Unlike the (local) nonlinear Schrödinger equation, which is integrable and can be solved via the inverse scattering transform method 1,3 , the NNLSE given by Eq. (1) is non-integrable 9,29 and cannot be solved analytically. In a non-integrable nonlinear system, chaos often appears. Chaos is generally agreed to denote the aperiodic long-term behavior of a bounded deterministic system that exhibits sensitive dependence on initial conditions. Scientific RepoRts | 7:41438 | DOI: 10.1038/srep41438 And the most common criterion for chaos is a positive Lyapunov exponent, which means that two initially arbitrarily close trajectories in phase space diverge exponentially in time [36][37][38] .
In this letter, we investigate the evolution of the (1 + 1)-dimensional NNLSE with the exponential-decay response function for the initial inputs of the HGT stationary solutions. As has been mentioned 20 , the HGT stationary solutions with the hump-number more than 4 always evolve unstably. We, however, find that such an unstable evolution of every HGT stationary solution can develop into a chaotic state, which is characterized by the positive Lyapunov exponent and spatial decoherence. Moreover, it also exhibits the soliton-like properties: the invariant statistic width during the evolution and the quasi-elastic collisions during the interaction. Therefore, we refer to these entities as chaoticons, as they are termed for the spatiotemporal chaotic localized states in the dissipative systems 39,40 . We believe it is the first time, to the best of our knowledge, to present the solutions in the conservative system (Hamiltonian system 36 ) described by the NNLSE which possess both the chaotic and soliton-like properties.

Unstable evolution of HGT stationary solutions
We consider here the NNLSE [Eq. (1)] with the exponential-decay response function 13,20,22 m m which has a singularity at x = 0. This case corresponds to the model for the propagation of the (1 + 1)-dimensional paraxial optical beam in nematic liquid crystals, which can be described by the coupled partial differential equations 10,11,15,16  where q is the dimensionless slowly-varying complex amplitude of the optical field, x and "time" t stand for, respectively, the dimensionless transverse coordinate and the dimensionless propagation direction coordinate, the "potential" − V represents the nonlinear induced refractive index. Obviously, the set of equations (3) is equivalent to Eqs. (1) and (2). The relative scale of the characteristic length of the response function w m to the statistic width of the wave w denotes the degree of nonlocality 13,14 , where w is defined by the second-order moment

N N
where u N is a real function and b N is a real constant. It was numerically found that in the case with the exponential-decay response function u N (x) is of N-humps (N = 1, 2, 3, … ) HGT-structure 20,41 , specially, u 1 (x) has a single-hump Gauss-type shape. It was also proved that u N (x)(N ≥ 2) can exist only when the parameters w m and b N satisfy . We simulate Eqs. (1) and (2) with the initial inputs of the HGT stationary solutions q(x, 0) = u N (x) by means of the split-step method 43 . The case of strongly nonlocal nonlinearity [w m = 10 and w(0) = 1 unless otherwise stated] is considered. The unstable evolution 41 of the HGT stationary solutions (N > 4)are given in Fig. 1, where only solutions with N = 7 and 12 are displayed without loss of generality. It is clear that the profiles starting from regular multi-humps turn to be irregular shapes, several of which are shown in Fig. 1(e). The evolution diagrams remind us the behavior of chaos.

Chaotic behavior: positive Lyapunov exponents
Since a positive Lyapunov exponent is a signature of chaos, we explore the maximal Lyapunov exponent [36][37][38][39] for the evolution of the HGT stationary solution. According to refs 44-47 the maximal Lyapunov exponent is computed by , which is the distance between two functions q 1 (x, t) and q 2 (x, t) in the Hilbert space (the L 2 norm in the Hilbert space), the two initial values q 1 (x, 0) = u N (x) and q 2 (x, 0) = u N (x) + r(x), and r(x) is a random perturbation function (as small as machine precision allows, e.g., in the order of 10 −8 ). For a partial differential equation, the Lyapunov spectrum (all the Lyapunov exponents sorted decreasingly) converges to a smooth curve, although the number of Lyapunov exponents is dependent on the number of discretization M 46 . In our numerical calculations, we choose M = 2048, 4096, 8192, each of them with the window size L = 40, 50, 60. We find that the statistical errors of the maximal Lyapunov exponents for every initial input are less than 10% (the maximal Lyapunov exponents are nearly independent on r(x) or the renormalization step-size). Then we get the average values of the maximal Lyapunov exponents for the HGT stationary solutions with N ≤ 12, as summarized in Fig. 2. We can see obviously that the maximal Lyapunov exponents for the unstable stationary solutions are all positive and increase monotonously with N, while those for solitons are equal to zero. The occurrence of chaos can be understood as a consequence of the complex interactions among humps. The more humps the profile possesses, the more complex the interactions are, which leads to a higher degree of the chaos.
It is especially important to make clear that the chaotic phenomenon described above is due to the intrinsic nature of the system but not numerically induced chaos 48 . Although the numerical method applied is not symplectic, it has been demonstrated that this is not the relevant issue for an infinite dimensional Hamiltonian system 45,49 . We have confirmed numerically that the scaling property of the maximal Lyapunov exponents match the transformation invariance of Eq. (1) 17 . The maximal Lyapunov exponents for the HGT stationary solutions with a given number of humps under the condition of = ∼ w w k / (k is a positive constant) and = ∼ w w k / m m satisfy λ λ =  k 2 within the error range allowed. The satisfaction of the scaling property is a stringent test for the reliability of numerical computations 45 .
Next, we will prove that the maximal Lyapunov exponents coincide with the growth rates of the initial numerical errors. Although, literally, the maximal Lyapunov exponent measures the typical exponential rate of growth of an infinitesimal perturbation, the growth of a noninfinitesimal deviation is usually well described in this way. The numerical error of the HGT stationary solutions computed by the Newton iteration method in double precision is assumed to be of the order of 10 −9 . It will make sense that the regular profiles of the initial HGT stationary solutions will be considered to become completely irregular once the deviation reaches the order of 1. We can, therefore, estimate t c (the critical time of becoming completely irregular) by = λ − e 10 1 t 9 c , thus obtain t c ≈ 20.7/λ, as shown in Fig. 3(a). From the other aspect, the process of turning to be irregular for the profiles can also be revealed directly in the evolution. Let's consider the skewness (or the third-order central-moment) of the intensity Obviously, there is s(0) = 0, since |u N (x)| 2 is symmetric. Figure 3(b) shows the evolution of s for the HGT stationary solutions with N = 7 and 12 in the time interval [0, 100] and [0, 50], respectively. It can be seen that |s| starts from the close neighbour of zero and then rises abruptly around a certain t cs , which is defined as . To a great degree, the boom of the skewness indicates the complete irregularity of the intensity profiles. That is to say, t cs represents the critical time of becoming completely irregular attained from the direct statistic method. The comparison between t c and t cs is given in Fig. 3(a). It is evident that the two curves always stay close to each other, which suggests that the critical times evaluated from the above two approaches agree approximately. Then we are certain that the maximal Lyapunov exponents obtained indeed indicate the exponential growth rates of perturbation.
In addition, it is expected that the spatial patterns will be spatially decorrelated in a system described by a partial-differential evolution equation with temporally chaotic behavior [50][51][52][53] . Then we calculate the spatial cross correlation function of two long enough wave-amplitude series at locations ξ and η  where the superscript * denotes the conjugate complex. The modulus of c for the HGT stationary solutions are depicted in Fig. 4, from which we can see that |c| equals 1 along the line ξ = η and decreases rapidly with the separation of two locations. The quick drop of correlation in the x direction means the spatial decoherence 52,53 .

Soliton-like property I: invariant statistic width
Chaotic as they are, the evolution of the unstable HGT stationary solutions maintain almost invariant statistic width w, as shown in Fig. 1(f). The standard deviation of w during the evolution of every HGT stationary solution is less than 0.02. It is well-known that 1,7,9 one of two intrinsic properties for the soliton is its invariant diameter  (width), thus we can conclude that the dynamic evolutions of the unstable HGT stationary solutions with invariant statistic widths are of such a soliton-property from the statistic point of view, even though their profiles during the evolutions are not constant. We can also see, as discussed next immediately, that they still possess the other soliton-property: a particle-like interaction. Because there co-exist the chaotic property and the soliton-like property during the dynamic evolution of the unstable HGT stationary solutions, we refer to them as chaoticons. Although we are not first to use the term "chaoticon", the intension in both mathematics and physics of the chaoticon here is completely different from that for the spatiotemporal chaotic localized structures in dissipative systems 39,40 .

Soliton-like property II: interaction of quasi-elastic collisions
Amongst all soliton properties, the important fascinating one is the particle-like interaction 7,8,23 . In order to check whether the chaoticons have such a property, we explore the interaction of two chaoticons that are initially identical and paralleled, which is presented in Fig. 5. The initial separation between chaoticons is large enough (8 times larger than w) to prevent the overlap of waves, and for each case of different Ns, both of the initial chaoticons are q(x, t 0 ) for t 0 ≥ t c , which means that the inputs are completely irregular states. We can observe that the two chaoticons attract each other, and then combine and separate quasi-periodically, much like elastic collisions between two particles. In fact, they will eventually fuse together accompanied by small energy loss to radiation after a much longer evolution. Hence, strictly speaking, the interaction is quasi-elastic.

Some remarks
Firstly, it is worth underlining that the evolution of the unstable stationary solutions in generally or weakly nonlocal nonlinearity 13 is entirely different from those in strongly nonlocal nonlinearity discussed above. In relatively weak nonlocality, the unstable HGT stationary solution will break up and form a set of single-hump profiles by emitting remnants of their energy, which is an unbounded state since the radiation waves arrive to infinity 27 . We have also found that a time when the wave begin to break up increases exponentially with w m for every HGT stationary solution with a given w. It means that the HGT stationary solutions in stronger nonlocality will evolve longer before they break up. Therefore, the radiation waves are believed to be absent in a strongly nonlocal nonlinear case 27,28 . Secondly, although the system considered is the (1 + 1)-dimensional NNLSE with the exponential-decay response function, our work may be readily extended to systems with different response functions 21 or even higher dimensions 26 . Thirdly, the mechanism that supports the localized spatiotemporal chaos in our conservative system is still an open problem, while the existence of dissipative chaoticons is explained analytically as the pinning and interaction of the fronts 40 .

Conclusions
We have found that the unstable evolution of the (1 + 1)-dimensional NNLSE with the exponential-decay response function for the initial inputs of the HGT stationary solutions will evolve into a new kind of chaoticon, which occur only in the case of strongly nonlocal nonlinearity. The chaoticon exhibits both chaotic and soliton-like properties. The chaotic behavior is signified by the positive maximal Lyapunov exponents and spatial decoherence. The soliton-like property is demonstrated by the invariant statistic width during the evolution, as well as the quasi-elastic collisions during the interaction.