A new kind of chaotic diffusion: anti-persistent random walks of explosive dissipative solitons

The solitons that exist in nonlinear dissipative media have properties very different from the ones that exist in conservative media and are modeled by the nonlinear Schrödinger equation. One of the surprising behaviors of dissipative solitons is the occurrence of explosions: sudden transient enlargements of a soliton, which as a result induce spatial shifts. In this work using the complex Ginzburg–Landau equation in one dimension, we address the long-time statistics of these apparently random shifts. We show that the motion of a soliton can be described as an anti-persistent random walk with a corresponding oscillatory decay of the velocity correlation function. We derive two simple statistical models, one in discrete and one in continuous time, which explain the observed behavior. Our statistical analysis benchmarks a future microscopic theory of the origin of this new kind of chaotic diffusion.


Introduction
Localized structures in dissipative systems are persistent solutions of homogeneous nonlinear wave equations with a spatial support that remains bounded as a result of a balance between dissipation and nonlinearity. These structures can show particle-like behavior similar to the solitons of conservative models. Thus, the name 'dissipative solitons' refers to localized structures with near-solitonic behavior [1][2][3].
One of the peculiar behaviors of dissipative solitons is the occurrence of explosions, which have been the subject of intense research in the last 20 years, as documented for instance in [4][5][6]. Although originally discovered in models of Nonlinear Optics such as the nonlinear Schrödinger equation (NLS) with dissipation and quintic nonlinearities, they have been found as well in other models in wide regions of parameter space. They can also appear in reaction-diffusion systems of multiple species if strong cross-diffusion is present [7]. More interestingly, explosions of solitons have been observed in several experiments of propagation of strong light pulses: in Ti:sapphire mode-locked lasers [8] and in double-clad ytterbium-doped fibers [9][10][11][12][13].
Generally speaking, explosions can be considered as the outgrowth of a fluctuating tail of a soliton that merges with the original soliton and becomes transiently a wide and large oscillating structure which then decays to a soliton that looks like the original modulo a spatial shift. Despite the similarity of the pre-explosion and post-explosion states, the lengths of the time intervals between explosions are irregular and no two explosions are perfectly identical.
It has been shown that explosions are the result of chaotic behavior, and its onset (as let us say a certain control parameter is varied) shares several features with the Ruelle-Takens route [14], torus-doubling bifurcation [15], and attractor-merging crisis [16]. Although explosions are the result of purely deterministic chaotic dynamics, the presence of additive noise can trigger explosions when they are not expected deterministically and make them more frequent when they do exist [17].
Here in this work, the focus will be the deterministic spatial diffusion of a dissipative soliton based on a generic model of nonlinear extended media that features an ordinary Laplacian for spatial coupling. In particular, we address the long-time behavior of solitons and the persistence of their zig-zag-like motion. A detailed study of the statistics of the soliton motion leads to a description of this kind of chaotic diffusion in infinite-dimensional systems as an anti-persistent random walk [18,19] thus bringing together the research fields of nonlinear spatially extended systems governed by partial differential equations and statistical physics of anti-persistent random walks. An extension of this random walk to a time-continuous model leads to further physical insight such as a deeper understanding of the velocity correlation function.

Explosions of dissipative solitons
In this article, we will use the cubic-quintic Complex Ginzburg-Landau (CGL) equation for modeling the dynamics of the envelope of the electric field, the complex amplitude A(z, t): This equation has been proposed in Nonlinear Optics as a model of passively mode-locked lasers [4][5][6][20][21][22]. Furthermore, the CGL is relevant in diverse fields such as the convection of binary fluid mixtures [23], chemical reactions on surfaces [24], and the Faraday instability in granular matter [25] and colloidal suspensions [26]. The left-hand side has the terms of the NLS, and the right-hand side includes a dispersion term with a second temporal derivative and nonlinearities associated to the Kerr effect. After interchanging space and time and some other trivial changes, we arrive at the more traditional form of the CGL that we will use hereinafter to facilitate the discussion of diffusive motion: and the constant parameters  m Î and , ,  d b g Î . This model is physically relevant for μ<0 and close to its subcritical bifurcation: when it has one stable zero solution and two harmonic non-trivial solutions, one unstable and another stable. The model has the following symmetries: translation in space, reflection in space, phase rotation, and phase conjugation.
The one-dimensional CGL, equation (2), was integrated from a localized initial condition A(x, 0) using a split-step Fourier method. The size of the domain L was chosen large enough so that the amplitude was truly localized: infinitesimally small outside the core of the soliton. To accelerate the simulations, we used graphical processing units programmed with the help of software PyOpenCL [27].
Depending on the specific values of the parameters, the pulse can show a variety of behaviors. Apart from being explosive, solitons can be stationary (in absolute value), radiating, and weakly chaotic. Here, we focus on a regime where explosions are present, μ c =−0.213 75<μ<0 (using the values given in figure 1 for the other parameters). The deterministic time evolution of an explosive solution of equation (2) is depicted in figure 1. As shown in the figure, there is a period of time when the soliton radiates small waves to both sides. Then, almost unexpectedly, one (either left or right) of the emerging waves grows, overcomes the original soliton, becomes wide and tall, then decays to a soliton that has roughly the original shape but is shifted in space. The process takes around 5 units of time, and the intervals between explosions are in the range between 20 and 50 units of time. The size of the spatial shift is comparable to the width of the soliton.
An explosive soliton is a stable solution, as opposed to other structures that may vanish or annihilate. For instance, for the same model, equation (2), holes were found to annihilate with each other, but during the transients, a rich anomalous dynamics appears [28]; normal diffusion of solitons in a Bose-Einstein condensate induced by heterogeneities has been reported [29]; transient anomalous diffusion of solitons has been found [30] in nonlocal random media. Figure 2 shows the evolution of a soliton for a longer time window including 6 explosions. There are some key observations from figure 2 that will guide the rest of the article: the times between explosions are not constant but have fluctuations; the directions of the spatial jumps seem to follow an almost perfectly alternating sequence; and the size of the jumps show a small variability. Interestingly, the times between explosions are such that the sum of two consecutive intervals is roughly constant.
Our observations and the perfect zig-zag patterns of localized structures [31] are qualitatively different. The one-dimensional and the two-dimensional behaviors of the soliton (the latter was reported in [32,33]) require separate treatments, since in two dimensions, the soliton can suffer perfectly symmetric explosions that do not induce spatial shifts.
In the present work, the main topic will be the long-time statistics of explosions and the diffusive motion of solitons. Therefore, the interest will be the overall motion induced by explosions and not the details of individual explosions or the times between consecutive explosions.  3. Statistical characterization of the soliton motion 3.1. Trajectories of solitons Despite its explosive behavior, a dissipative soliton remains a localized structure with a well-defined center for unlimited periods of time in the absence of external perturbations. Thus, the question of a possibly diffusive motion of the center's trajectory is perfectly meaningful. As explained in the previous section, the intrinsic chaotic behavior of a soliton is manifested in an apparently random sequence of jumps.
A common observation was that when starting from a perfectly symmetric initial amplitude A(x, 0), the first few explosions were symmetric and looked as if both tails grew at the same time. But this situation is transient and highly non-generic, in fact, eventually, all subsequent explosions are asymmetric inducing either a displacement to the left or to the right.
A typical intensity pattern as function of time is shown in figure 2. To characterize the motion of a soliton, we define its 'center of mass' in a periodic domain as This will allows us to generate trajectories of solitons starting from localized roughly Gaussian initial conditions centered at x(0)=0.  Figure 3(a) shows a sequence of 8 perfectly alternating explosions. Despite the regularity, small fluctuations in the inter-explosion times δτ k and the jumps δx k are visible. Figure 3(b) shows a defect in a long zig-zag pattern. These defects appear unexpectedly and cannot be predicted from a localized initial condition. In figure 3(c), we show the deterministic evolution of ten slightly different initial conditions A(x, 0) for μ=−0.180. The total integration time was T=10 6 and each realization included roughly 40 000 explosions organized in long zig-zag patterns that result in small net displacements. After a few explosions, the centers x(t) describe trajectories that look like random walks. Comparing trajectories obtained for the same parameter values, we observe that all trajectories are different but share some common qualitative features.

Squared displacements 3.2.1. Ensemble-averaged squared displacements (EASDs)
The basic characterization of diffusion is given by the asymptotic time dependence of the mean-squared displacement, which can be determined as an ensemble average or as a time average. In the former case, squared displacements are averaged over an ensemble of independent realizations: resulting in the EASD. Figure 4 shows the EASD obtained for 800 independent realizations and a long simulation time. At first sight, the EASD depends linearly on time: (12)). Interestingly, on a finer time scale, the EASD shows oscillations (displayed in the left inset in figure 4) which vanish for larger values of τ in the limit of N  ¥ (see the right inset in figure 4). As we will show, these oscillations are a result of the alternating sequence of spatial shifts.

Time-averaged squared displacements (TASDs)
A complementary perspective of diffusion can be obtained by analyzing individual trajectories and their temporal structure using the TASD for a lag time τ<T: In figure 5, we show the corresponding TASD for each trajectory. It shows two interesting features. First, a large scatter is observed that grows as larger lag times τ are considered. Second, each curve shows oscillations that are shown in the inset. We found similar behaviors for other values of μ in the range μ c <μ<0.
The large scatter of the TASDs is an effect of finite simulation time T and is observed even for Brownian motion. As we will show later, the scatter of the TASD is actually consistent with Brownian motion. As longer and longer simulations are used, the variance of Despite the wild scatter between trajectories, the ensemble average of the TASDs (EATASD): does grow linearly with τ parallel to the EASD as can be seen from figure 6. From these figures, we can estimate a diffusion coefficient of the soliton, D 0.000 57 The offsets of the EASD with initial aging time and the EATASD roughly coincide and can be explained in the context of our model in section 4.1. The offset of the EASD without aging is larger due to the transient behavior of the soliton motion in the beginning.
As already indicated, large trajectory-to-trajectory variations lead to (at least apparent) 'irreproducibility' of statistics of individual trajectories despite long measurement times T. As x 2 T t áD ñ ( ) takes wildly different values In order to answer this question, we analyze the distribution p(ξ, τ) of normalized TASDs: In figure 7(a), we show a histogram for the random variable x t ( ) obtained from individual trajectories, and we compare it with the corresponding analytical result for Brownian motion. According to Andreanov and Grebenkov [34], the latter can be approximated by a generalized Gamma distribution: with three parameters a0, b>0 and  n Î , which depend on the ratio τ/T, and K ν denotes the modified Bessel function of the second kind. Assuming that the ratio τ/T is small, Andreanov and Grebenkov obtained for the parameters: We á ñ = ´and the intercept is c=5.445±7.938. This estimate was obtained by fitting the full data range 0τ10 6 . Eliminating the first 1000 units of time, the estimate for the slope is D 0.000 5617 3.06 10 1 5 á ñ = ´and for the offset c=0.519±7.584. The uncertainties were computed using the bootstrap method. The left inset shows a narrow time window of the same curve with oscillations that are a result of the antipersistence of spatial shifts. As it can be seen in the right inset, the amplitude of these oscillations (measured as the square root of the mean square deviation with respect to the linear fit) scales with N 1 (indicated with a line of slope 1 2 in the logarithmic plot) and vanishes for larger values of τ in the limit of N  ¥. used this approximation for the theoretical curve in figure 7(a), which is in good agreement with the distribution p(ξ, τ) for the diffusing solitons.
The variance of x t ( ) is known as the ergodicity breaking parameter. For Brownian motion, it only depends on the ratio r=τ/T and is given by [34]    In figure 7(b), we compare this theoretical result with the ergodicity breaking parameter for the diffusion of the solitons. Obviously, if the ratio τ/T grows, the variance of x t ( ) also grows. For a fixed value of τ and an increasing value of T, the distribution p(ξ, τ) becomes narrower and concentrated at ξ=1, i.e. the TASD of a single trajectory will agree with the EASD.

Distribution of generalized diffusivities (DOGD)
Another tool that can be used to analyze diffusion processes is the DOGD [33,[35][36][37][38][39]. A generalized diffusivity is thereby defined as a single squared displacement which is rescaled with the asymptotic time dependence of the mean-squared displacement: D . Analogous to the mean-squared displacement, in order to obtain a distribution, the required squared displacements can be obtained from an ensemble of diffusing particles or from one long single-particle trajectory. If the scaling exponent α is equal to the exponent of the asymptotic increase of the corresponding MSD, the mean generalized diffusivity is equal to the generalized diffusion coefficient. Therefore, the DOGD captures the fluctuations during the diffusion process around the generalized diffusion coefficient obtained from the mean-squared displacement. The choice of the exponent α and the time frame τ may reveal non-trivial pieces of the process dynamics.
The DOGD obtained as ensemble average is defined by and captures fluctuations over the various members of an ensemble. The DOGD obtained as time average of a single realization is defined by and captures fluctuations along a single-particle trajectory. This distribution is random due to the random nature of its mean D T t á ñ a ( ) . The distribution varies from one trajectory to another. We consider its ensemble average p D, ) in order to obtain a non-random distribution, which can be compared with analytical results.
For Brownian motion, and using α = 1, the distributions p 1 E (D, τ) and p D, ) coincide due to the ergodicity of normal diffusion with respect to the squared displacements and are given by the formula: which is a rescaled version of the propagator and does not depend on τ due to the correct choice α=1. The first moment of p 1 (D) is In figure 8(a), (b) one can see a good agreement (in semi-logarithmic plot) of the analytical expression p 1 (D) and the histograms extracted from numerical simulations.

Fine-scale motion of solitons
Sampling the trajectory at the midpoint between two consecutive explosions, one can extract a series of shifts δx k and a series of time intervals δτ k . Figure 3(a) shows the meaning of δτ k and δx k . Although these two series include a lot of regularity, the presence of fluctuations is key to understand the long-term diffusive motion of the soliton.

Spatial fluctuations
Although the shifts show a zig-zag pattern for very long intervals, there are two imperfections that contribute to the diffusion process. First, there are some defects of the zig-zag pattern where two or more successive jumps follow the same direction. Second, the sizes of the jumps are not always equal. Figure 9 shows the distribution of the values that δx k takes. The possibility of left and right jumps and the fluctuations in jump sizes lead to two peaks with finite widths. Although small, these fluctuations contribute to the long-term diffusive motion of the soliton. Figure 10 shows the absolute value of the correlation function of the spatial shifts and of the algebraic signs of the spatial shifts. Both curves decrease exponentially with the same slope in semi-logarithmic plot indicating that the correlation of the spatial shifts is only due to the correlation of the algebraic signs and not due to a correlation of the sizes of the spatial shifts. This is confirmed by figure 11, which indicates a correlation length of the fluctuations in the size x k k d = ℓ | | of the spatial shifts smaller than one. Therefore, the size of the spatial shifts can be considered as uncorrelated.  The correlation of the algebraic signs is due to the zig-zag pattern of the spatial shifts and the occasional defects. Assuming that the direction of the kth jump only depends on the direction of the (k−1)th jump, we can define a persistence parameter (probability) a that two successive jumps go in the same direction. Accordingly, the probability of a zig-zag pattern of length Δk without any defect is The exponential decrease is consistent with figure 12. We numerically estimated the persistence parameter, a=0.000 86, which is roughly the same as the exponential decay in figure 12, 0.000 76. Note that for small values of Δk, the numerically determined probability is larger than expected. This can be explained by the fact that we observed that a defect (two jumps in the same direction) is often immediately followed by another defect (for instance, three consecutive jumps in the same direction). This is the reason why the mean value of Δk is a little smaller than expected from the exponential decay. ) explosions for a=0.000 86 is smaller than the numerically determined value in figure 10, 870 explosions. A possible reason is that a certain number of defects occurring as several pairs leads to a longer correlation length as if those defects would be equally distributed over the trajectory. However, we want to emphasize that the simple assumption of the persistence parameter a leads to theoretical models, discussed in section 4, which can excellently reproduce the diffusive properties of the soliton motion such as the diffusion coefficient D 1 á ñ.  ) indicates that the correlation length is equal to 0.69<1, so this correlation can be safely neglected in the randomwalk models presented in this work.

Fluctuations in the timing of events
The distribution of the inter-explosion times δτ k by itself, as long as it has a finite mean value, does not essentially influence the diffusion process. Only the mean value of the distribution determines the diffusion coefficient. Figure 13 shows a numerically determined distribution which is centered around a well-defined mean. Figure 14 shows the fluctuations of the inter-explosion times δτ k around the mean 26.27 dt á ñ = . Looking at a smaller time scale, it is possible to observe a zig-zag pattern with a slow modulation. The short-long-short-etc sequence was already identified in figure 2. The figure suggests 52.54 There is a slow modulation of the inter-explosion times: δτ 2k takes values above the average for a few explosions, then below the average for some others. Every time the previous curve changes sides shows a breaking in the pattern shortlong-short-etc. Figure 15 shows the autocorrelation function of the inter-explosion times δτ k , which is oscillating due to the short-long-short-etc pattern with an exponentially decreasing envelope.

Velocity correlations
An illustrative way to analyze the anti-persistence of the spatial jumps is the use of the velocity autocorrelation: (first averaged over individual trajectories and then over the whole ensemble) where v t x t = ( )˙( ) is the velocity of the center of mass of the soliton. Figure 16 shows this correlation function consisting of a sequence of sharp positive and rounder negative peaks, both decreasing in amplitude for increasing τ. Not surprisingly, the velocity autocorrelation shows the same persistent oscillations as the EASD, because they are connected by the Green-Kubo relation:      (6)) T=100 0000 and then over the N=800 realizations. The rounder negative peaks are induced by the anti-persistent spatial shifts and the fluctuating inter-explosion times (between 20 and 30 units of time); in fact it is the sum of two consecutive inter-explosion times, the quantity that is maintained almost constant (around 52.54), that explains why the first positive peaks are so sharp. The inset shows (in logarithmic scale) the decay of the negative and positive peaks. The slope of the exponential decay (depicted in red) is s=−3.97×10 −4 (estimated for 1000τ3000) that indicates that correlations last for around s 96 that can also be written as and approaches the diffusion coefficient D 1 á ñ in the limit of large τ as illustrated in figure 17.
There are several observations from the velocity autocorrelation function. General oscillatory behavior is a result of the alternating nature of x(t). Adjacent peaks with same algebraic sign separated by 2 52.54 t dt D = á ñ » correspond to mean separation of explosions with equal directions. Adjacent peaks with different algebraic sign separated by 26.27 t dt D = á ñ » correspond to mean separation of explosions with opposite directions. As time lag τ grows, two different effects take place: sharp positive peaks, a result of the almost constant sum of two consecutive inter-explosion times, become broader after around 20 explosions, and the overall autocorrelation function becomes more and more triangular. The correlation length of the velocity is roughly 96 explosions, a value that is much smaller than the correlation length of the spatial shifts in figure 10. The reason is that contrary to the sequence of spatial shifts, the velocity correlation is also influenced by the fluctuations of the inter-explosion times. This is confirmed by numerical simulations regarding our continuous-time model in section 4.2 (see also figure 19). The fast decay of the inter-explosion times' correlation is observed as a progressive broadening of the positive peaks in the velocity autocorrelation function.

Statistical models of soliton motion 4.1. Discrete-time models
The soliton motion x(t) can be modeled as where θ(t) is the Heaviside step function and t k is the time of the kth explosion: Because the diffusion coefficient in only influenced by the mean value of the inter-explosion times δτ k , the time between explosions can be set to this constant value, and, therefore, we can consider a time-discrete model: Now we separate the increments δx k into their directions s k =±1 and their sizes x 0 | | , so that the discrete equation of motion becomes The sequences s k { } and k ℓ { } are in general realizations of two stochastic processes, which may not be independent and may also be correlated in time k. We can now integrate the equation of motion to obtain  To proceed, we have to make assumptions about the two stochastic processes. The first assumption is the independence of the direction and the step size process so that where we can identify a diffusion coefficient, an offset, and an oscillating term. The diffusion coefficient is given by We want to emphasize that we investigated the soliton dynamics in the whole range of the parameter μ where explosions of dissipative solitons occur. We used steps of Δμ=0.01 and found that the anti-persistent random walk picture is correct, i.e. the diffusion coefficient is determined by equation insofar as one starts to observe a kind of superposition of two anti-persistent random walks, which also emphasizes the importance of the understanding we present in this article.

Continuous-time models
Although the fluctuations of the inter-explosion times δτ k around the mean dt á ñ, summarized in figure 14 and described in the previous section, cannot change the overall diffusion coefficient D 1 á ñ of the random walk, there are fine details that can be captured by including these fluctuations. Especially, they can explain the discrepancy between the correlation length of the velocity and the correlation length of the spatial shifts. Therefore, in the following, we incorporate the fluctuating inter-explosion times in our discrete-time random walk model and obtain, therefore, a continuous-time model. The starting point of the new model is the observation of the anti-persistence in the sequence of interexplosion times, presented clearly in figure 15. A sequence of waiting times between jumps that follows the correlation function in figure 15 can be numerically generated by a slightly modified AR(1) process: where ξ k is Gaussian white noise with 0, 1 , . 33 , the mean value of this process is asymptotically given by By choosing appropriate values of the parameters m, b, and σ, this slightly modified AR(1) process is able to reproduce the statistics of the observed inter-explosions times for the soliton motion. For instance, the autocorrelation function, equation (35), represented in figure 18 captures the detailed features present in figure 15.  The time-continuous model reproduces the velocity correlation function of the soliton motion, in particular, the progressive broadening of the positive peaks and the effect that the correlation length of the velocity is much smaller than the correlation length of the spatial shifts as depicted in figure 19.
Finally, the time-continuous model reproduces the EASD and the ensemble average of the TASD. Figure 20 shows both the EASD and the EATASD that are consistent with figure 6. The fit gives now D 0.000 584 1 á ñ = with an offset 2.8. The EASD (blue curve) even shows the persistent oscillations associated to the finite number of samples. The dependence of the EASD on the aging time t a , responsible for part of the vertical offset in figure 6, does not appear in the model, because it is due to the transient behavior of the soliton dynamics, which is not captured by the model.

Conclusions
Explosions of dissipative solitons of the CGL were already understood as an outcome of their chaotic behavior. Here, we have addressed for the first time the random and strongly anti-persistent behavior of the center of mass of the solitons and the particular features of the corresponding diffusive motion. Additionally, we have studied the anti-persistence of the times between explosions and its effects on second-order statistics such as the velocity autocorrelation.
There are several ways to analyze diffusive behavior in a deterministic context. Here, we integrated an ensemble of independent realizations of the soliton initialized with small variations. Then, we used the EASDs and the TASDs to characterize the time evolution of an ensemble of trajectories and the temporal structure of the trajectories. We obtained that the TASDs of the individual realizations have a wide distribution, even if their ensemble average for long time lags converges to the EASD that grows linearly with time.
The large variation among the TASDs of the individual realizations is an effect of finite simulation time T that is even observed for Brownian motion and gets reduced as longer trajectories are considered.
We presented two simple statistical models for the motion of the solitons. First, we presented a time-discrete random walk model which neglects the fluctuations of the waiting times but incorporates the fluctuations of the jump sizes and the defects of the zig-zag patterns. Both mechanisms are responsible for the overall diffusive properties, captured by the temporal behavior of the ensemble-averaged and the TASDs.
Second, we presented a time-continuous model which also incorporates the fluctuations of the interexplosion times. This model shows all the diffusive properties of the first model and reproduces further features of the stochastic process such as the velocity correlation function.
In summary, we have shown that a prototypical spatially extended system can show deterministic chaotic diffusion of localized solutions, which can be very well characterized as an anti-persistent random walk. We performed an exemplary thorough statistical analysis, which sets benchmarks for a future microscopic theory of the origin of this new kind of chaotic diffusion.