Multi-delay complexity collapse

Increasing the number of delays in nonlinear dynamical systems is generally assumed to lead to higher complexity, but “distributed delay” systems with an inﬁnite number of delays to lesser complexity. This paradox is studied here using the Lang-Kobayashi laser model by ﬁrst extending a recent method for Lyapunov exponent estimation from single to multiple delays. The Kolmogorov-Sinai entropy, permutation entropy, and time delay signature suppression initially increase with number of delays as the dynamics become more hyperchaotic. At a large number of delays that depends on the feedback strength, this trend reverses, leading to simpler dynamics. The phenomenon is also found in other delay equations, such as the Mackey-Glass system. A similar collapse is uncovered in the distributed delay case with broadening delay kernel.

Delay Differential Equations (DDEs) are a class of dynamical systems governed by some memory of the past.They exhibit a rich array of behaviors in a broad range of applications [1][2][3].Their multiple timescales [4] can be exploited for the generation of oscillatory [5] and chaotic time series, including the recently discovered laminar chaos [6], as well as for stabilizing complex dynamics [7] and unstable steady states [8].What role plays the number of delays M? Adding a second delay has been shown experimentally to induce hyperchaotic behavior [9].More recently, fiber gratings in semiconductor lasers that synthesize large numbers of random delays have been used to produce strongly hyperchaotic dynamics with application to random number generation and secure chaos communication [10].
In contrast, stabilizing effects can also arise with multiple delays.Second and third delays added to the Lang-Kobayashi (LK) laser system can produce continuous-wave dynamics with enhanced stability [11].A second delay can stabilize unstable periodic orbits in the Mackey-Glass (MG) system [12], and many feedback terms with delays that are integer multiples of a fixed delay can control unstable periodic orbits [13].Feedback control of Chua's circuit with two different delays can replace chaos by simpler stabilized steady states [14].Dependence on an infinite number of past states defines distributed delay systems with integrodifferential dynamics that are also generally "simpler"-i.e., with a less positive spectrum of Lyapunov exponents (LEs) and smaller dynamical entropy especially when the feedback kernel is broad (see, e.g., Ref. [15]).
Although details matter to some extent, a general paradox does emerge: More delays can lead to stronger chaos, but a large number of delays-and sometimes even two delayscan lead to simpler dynamics.This M-dependent transition is not well understood, and our paper studies this paradox using a semiconductor laser model with many delays.Such lasers with optical feedback are useful in our context due to their ability to produce high-dimensional hyperchaotic behavior.Modes with frequency less than the characteristic frequency of the laser in isolation (i.e., without feedback) can destabilize the fixed point through a Hopf bifurcation, yielding a limit cycle, then a torus, and eventually chaos and hyperchaos [16].
Information about the presence of a time delay can be extracted from the chaotic time series by the autocorrelation function (ACF), delayed mutual information or permutation entropy [17].Setups involving optical injection, rings [18], or phase conjugate optical feedback mirrors [19] have been shown numerically and experimentally to suppress the time-delay signature (TDS) in the ACF, thereby increasing unpredictability.TDS suppression can be achieved with multiple random (but fixed) delays, a useful property for random number generation [10,20].To understand the complexity transition at higher M, we make use of the TDS from the ACF, but also estimate dynamical invariants, including Lyapunov exponents through an extension of a method for a single delay [21].The dimensionless laser Lang-Kobayashi semiconductor laser model with multiple discrete delays can be expressed as where E (t ) is the slowly varying complex electric field, N (t ) is the carrier density, M is the number of delays, and κ is the feedback strength.The τ i are the delay times.The intrinsic relaxation oscillation period of this system is τ ro = 0.75.The values and descriptions of the parameters used in Eq. ( 1) are presented in Table I.We assume that the random fiber grating consists of M segments, each with length half of L = 1.8 mm.The delays are picked randomly for each segment as where i is from one up to the maximum number of delays M, with τ avg > (M + 1)L/(2 c) is randomly chosen from [0, L c ], the average delay is τ avg = 5 and c = c T where c is the speed of light and T normalization constant for time which is 1 ns.Time series were obtained using a fourth-order Runge-Kutta routine with a time-step 4 × 10 −4 .Figure 1 shows the interplay of the different timescales with M for fixed feedback strength κ.For the single delay case, the ACF displays a sharp peak at multiples of the delay [Fig.1(b)].For 20 delays, Fig. 1(e) displays no time-delay signature around the time delays, however, peaks around τ ro become more pronounced.
For large M, such as 160 delays in Fig. 1(g), the attractor appears like a torus, and the ACF displays two components with timescales of τ avg and τ ro .This is reminiscent of the single-delay weak feedback case where the external cavity causes periodic oscillations.The time delays still affect the dynamics; their ACF signatures are seen around lag τ avg with contributions of similar magnitude to that of the intrinsic response.M = 180 leads to a limit cycle with period τ ro , and only the fast intrinsic frequency is left.Increasing M to 1000 leads to a fixed point in Fig. 1(m).
The central peak at lag 0 broadens as the number of delays increases, and it is modulated by oscillations with a period near the intrinsic period τ ro (Fig. 2).It is known [22] that such oscillations are visible in ACF peaks when the feedback strength is not too high; otherwise the ACF peaks are too narrow to resolve those oscillations.In this sense, increasing the number of delays to 20 seems to have a similar effect to lowering the coupling strength.But, as we will see below, the complexity measures behave nonmonotonically as a function of M in this range of coupling strengths.
We can rewrite Eq. (1a) as a function of optical phase φ and intensity FIG. 3. Probability density of the total feedback term in Eq. (3a) for different number of delays.
The numerically estimated densities of the feedback term in Eq. (3a) are shown in Fig. 3.The standard deviation clearly decreases with increasing M, and the shape goes from unimodal to multimodal.The density appears almost skewed Gaussian for M = 20, but this is generally not the case.
A small value of the ACF at a lag around the mean time delay is suggestive of a parameter regime that is more strongly chaotic and vice versa.This intuition arises because this simple ACF measure is inversely proportional to the TDS, which generally decreases the stronger the chaos is-although no direct theoretical link has been established between low TDS as measured in the ACF and positive sub-Lyapunov exponents (see Ref. [23]).We thus plotted the maximal peak of the ACF around the mean delay for different M and κ values in Fig. 4(a).Blue (dark gray) regions, starting at M = 1 and κ = 3, correspond to regimes where TDS is small, whereas red (intermediate gray) regions, starting at M = 5 and κ = 2, reflect torus, limit cycle, and fixed point behavior.Increasing κ at fixed M, first the value of this maximal peak decreases, corresponding to transitions from fixed point to limit cycle to chaos, until it reaches a minimum.In this low-signature regime, the dominant time scale is the intrinsic response of the laser.For larger M, however, a stronger feedback coefficient is needed to enter the chaotic regime, suggestive of relatively simpler dynamics for large M.
When the number of delays become infinite, the state of the system depends on a continuum of past states.Such ecological [24] or neural feedback [15] distributed delay differential systems are known to have simpler dynamics, especially when the delay kernel is broad.They can be expressed using integrodelay differential equations, where g(τ ) denotes the kernel which must satisfy Here, we used the uniform distribution for g(τ ), F [N (t ), E (t ), E (t − τ )] is the right-hand side of Eq. (1a), and N (t ) evolves as in Eq. (1b).The maximum peak of the ACF near τ avg for different widths of the delay distribution σ and feedback strength are depicted in Fig. 4(b).The periodic structures in Fig. 4(b) show that semiconductor laser dynamics with distributed delays keep transitioning between states of high and low TDS, until no more chaos exists (σ 0.12).These structures occur when the width of the distribution (σ ) times ω is approximately 2nπ where n is a positive integer.In this case, the coupling between delayed terms and the current state requires a very strong feedback coefficient since the sum of the exp(ωτ i ) factors is a very small value.Figure 5 shows the collapse of complexity for this distributed LK model for increasingly broad delay kernels, i.e., for a larger range of delays.Note that simulating these dynamics involves a fixed discretization of the delay kernel set here to the integration time step; the trapezoidal rule was used for the integral.
Calculation of the Kolmogorov-Sinai (KS) entropy can generally be used to quantify the unpredictability of the motion since it measures the average information loss rate.It can be computed by summation of all the positive Lyapunov exponents, but the estimation of the Lyapunov spectrum for multiple delays is an open problem.Here, we resolve this issue by converting the multidelay DDE to m t + 1 ordinary differential equations [21], where the u j (t )'s are functions defined on the Chebyshev nodes θ i j in the interval between x(t ) and x(t − τ M ), f [x(t ), . . ., x(t − τ M )] is the right-hand side of Eqs.(1a) and τ M is defined in Eq. ( 3).The quantity m t is defined further below.The m i + 1 Chebyshev nodes on each of the M intervals (τ i−1 , τ i ) between successive delay pairs in the discrete delay distribution are defined as This formula shows that the nodes between each pair of delays are spaced differently in a manner that depends on the m i values.It should be noted that only the evolution of u 0 (t ) depends on the delayed variable values x(t − τ i ).The d k j 's are the elements of the Chebyshev differentiation matrix.There are m i + 1 nodes between delays τ i and τ i−1 .The m i values can be determined by numerically solving [25] ) in such a way that M i=1 m i = m t and τ 0 = 0. We calculate the Lyapunov exponents from the evolution of the state vectors for the associated linearized differential equations of Eq. ( 6) with reorthogonalization after each time step using the modified Gram-Schmidt method [26,27].
Lyapunov exponents have not been computed previously for the multidelay LK model using this method, nor have they been estimated experimentally, and thus it is not possible to independently validate our new method for estimating Lyapunov exponents.We can, nevertheless, get a good degree of confidence by checking that it reproduces known exponent values when restricted to the single delay case.We plotted the four largest Lyapunov exponents in Fig. 6(a) for the LK model with a single delay in the weak chaos regime (at high feedback strength).We can see the dependency of the Lyapunov exponents on the time delay.As mentioned in Ref. [23], using a larger delay makes the Lyapunov exponents smaller such that after a limit, λ m τ becomes almost constant where λ m is the maximal Lyapunov exponent.The effect of feedback strength κ on these maximum Lyapunov exponents is shown in Fig. 6(b).The highest value of the maximal LE is found for κ = 13.Further increasing the feedback strength decreases these four maximal Lyapunov exponents, and they also get closer to each other.Parameters used in Eq. ( 1) for Figs.6 and 7 are taken from Ref. [28] in order to compare the values obtained for the LK model.We found a good match between the Lyapunov exponents calculated in Ref. [28] and we further extended the simulation for higher feedback coefficients.Namely, to compare some values, we got λ m = 0.44 and λ m = 2.43 for κ = 3.4 and κ = 7.8, respectively, whereas the values found in Ref. [28] were approximately 0.43 and 2.5.
We can obtain further validation using another test where a second feedback pathway with a different delay (τ 2 = 1.4) is added to the single delay case above (τ 1 = 2) and check the behavior of the exponents as the strength of the second pathway becomes small.Specifically, instead of normalizing the feedback strength in Eq. ( 1) by the usual constant factor 1/M, we assumed that each delayed feedback pathway has its own weight w i which is multiplied by κ.We performed the calculation for a two-delay system.In the first case, we set w 1 = w 2 = 0.5.We see, in Fig. 7(a), that the largest Lyapunov exponent occurs for κ = 19, that this value is larger than in the single delay case Fig. 6(b) and that the exponents again converge for larger values of κ.In the second case shown in  6(b).This is one validation of the Lyapunov exponent estimation method: the spectrum for two delays converges to that for the larger delay on its own when that larger delay pathway dominates.FIG. 8. KS entropy (a) as a function of the number of delays when the feedback strength is fixed, and (b) as a function of the feedback strength when the number of delays is fixed.Fig. 7(b), we weaken the second delayed feedback by setting w 1 = 0.95 and w 2 = 0.05.The maximum Lyapunov exponent now occurs for a smaller κ; in fact, the curve is seen to converge to the one in Fig. 6(b).In other words, when the shorter delayed feedback pathway makes a very weak contribution to the dynamics, the system behaves dynamically very closely to the single delayed LK model in Fig. 6, an intuition that is reflected in our exponent estimations.
The KS entropy is plotted versus M and κ in Fig. 8 using parameter values in Table I and τ avg = 5.One sees that for low κ, increasing M brings the KS entropy to zero as the dynamics transition from hyperchaos to limit cycle.There is a small remnant of chaos even for M = 16 for which the entropy is 0.3 (κ = 2); it is essentially zero for M = 40.For larger κ, the entropy is still relatively high at M = 40 but exhibits the key decreasing trend that qualitatively agrees with Fig. 1.Numerical calculation of entropies for the large M values in Fig. 1 are prohibitive, but the slow decay seen in Fig. 8(a) suggests that simpler dynamics may, in fact, consist of a small scale chaotic attractor with a simple (e.g., limit cycle) skeleton.For a fixed number of delays, increasing κ first brings on a strong chaos regime (unpredictability increases), followed by a weak chaos regime (unpredictability decreases) [Fig.8(b)].These notions were first introduced in Ref. [23]: for strong chaos, the nonlinearity due to the relaxation oscillation is high, and for weak chaos, the feedback strength is either high or low.It can be seen that for more delays, a higher KS entropy can be obtained.However, the starting entropy values for small κ shift rightward for larger M, i.e., predictability increases, and the motion is simpler.
We also calculated permutation entropy (PE) of the intensity time series to study the dependency of the complexity (entropy production) on M and κ [29].Since the mean time delay is approximately the same for any M, we used the time delay for the single delay case in Fig. 1 as the embedding delay.PE can take values between 0 and maximal complexity 1.In agreement with the analyses up to now, Fig. 9(a) shows that PE is a nonmonotonic function of M when κ is sufficiently large, although the decrease is very slow at large M. We can see in Fig. 9(b) that, the smaller M gets, the weaker κ must be to reach high complexity.But the ordering of the curves with M inverses at stronger feedback like for the KS entropy, with a stronger κ being needed to reach higher values of PE.Whereas increasing M in the discrete delay case or the kernel's width in the distributed delay case, the mean and standard deviation of the feedback term [last term in Eq. (3a)] decreases (Fig. 3).This is not surprising for the LK system where the sum of harmonic functions at random phases in the numerator is overtaken by the 1 M normalization factor.As M increases, the density of the feedback fluctuations that was unimodal at low M becomes more U-shaped, similar to that for a sine wave, and narrows as well.The intuitive picture offered by the complexity collapse is one where large M delayed dynamics increasingly resemble those of a periodically driven, nondelayed, and underdamped stochastic system with small residual chaotic fluctuations thickening apparently simple limit cycles [Fig.1(g)]; this analogy will be pursued elsewhere.
We also observed the complexity collapse with M in the Ikeda model after an initial increase (not shown).Next, to illustrate the generality of this multidelay complexity collapse phenomenon, we turn to its investigation in a nonlinear delaydifferential system where the flow does not contain harmonic functions, in contrast to, e.g., the LK and Ikeda systems.Specifically, we consider the MG model governed by Here, we set γ = 1, β = 2, and τ avg = 17.The discrete delays are chosen in the same way as for the LK model above: new delays required for a higher value of M are chosen randomly from the distribution and added to the delays used in the simulation for the previous lower value of M. We see, in Fig. 10, that, for a sufficiently large M, simple oscillatory behavior arises.The feedback function then has a multimodal stationary density, and the correlation time of the feedback function, i.e., the first time when the ACF of the feedback equals 1/e of its maximum at lag 0, increases.Figure 11  entropy that is compatible with what is seen in the phase plots of Fig. 10.In contrast to the LK model, however, the dynamical complexity in the MG system does not increase upon adding delays starting at M = 1 or, at least, not for these parameter settings.
In conclusion, we found a limiting number of discrete feedback delays for the LK model beyond which the complexity gradually collapses to simpler motion, characteristic of distributed delays with sufficiently broad memory kernels.Adding delays causes more interference pathways for the feedback light, enhancing the relaxation oscillation, and lowering the TDS.This oscillation eventually bifurcates to a stable fixed point.This phenomenon appears as a generic feature of multidelay dynamical systems and likely of networks with large numbers of delays, and its origin may be linked to the low effective dimension observed in stochastic distributed delay systems [30].Multidelay nonlinear systems may, thus, paradoxically exhibit simple behavior.A differential delay system owes its complexity, in part, to the continuous perturbations that the nondelayed system receives from delayed feedback at a specific time in the past relative to the present, i.e., to the underlying dynamics of a discrete-time nonlinear map.Adding more delays at different times amounts to attenuating this discrete nature of these perturbations thereby causing a reversion to simpler dynamics.Our paper highlights this progressive loss of complexity, which can, in fact, include steep drops as in the MG case.

FIG. 4 .
FIG. 4. (a) Impact of feedback strength κ and number of delays M on the maximal peak height of the ACF around the mean time delay for multiple discrete delays.(b) Impact of κ and width of the delay kernel on the maximal peak height of the ACF around the mean time delay in the distributed delay case.

FIG. 6 .
FIG. 6.The four largest LEs (a) versus the time-delay τ with fixed feedback strength κ = 35, and (b) versus κ for fixed time-delay τ = 2.The LK model here has only a single delay.

FIG. 7 .
FIG. 7. The four largest Lyapunov exponents versus feedback strength for a configuration with two delayed feedback pathways.The delays are τ 1 = 2 and τ 2 = 1.4.(a) Both pathways have the same weight: w 1 = w 2 = 0.5.(b) Pathway weights are asymmetric with w 1 = 0.95, w 2 = 0.05.The upper curve in (b) is close to the upper curve in Fig.6(b).This is one validation of the Lyapunov exponent estimation method: the spectrum for two delays converges to that for the larger delay on its own when that larger delay pathway dominates.

FIG. 9 .
FIG. 9. Permutation entropy (a) as a function of the number of delays when the feedback strength is fixed, and (b) as a function of the feedback strength when the number of delays is fixed.
(a) shows the behavior of the correlation time versus M.After M = 40, a monotone increase in the correlation time with M is observed until around M = 705 at which point it begins a rapid drop by about 10%. Figure 11(b) displays the behavior of the permutation entropy versus M, and a sharp drop is also seen around M = 705.The dynamics beyond the drop display a moderately shorter correlation time and a significantly lower

FIG. 10 .
FIG.10.Attractor evolution of the Mackey-Glass system for an increasing number of delays.Here, x(t ) describes the current state, whereas x(t − τ avg ) refers to the state at the average time delay in the past.Simplification of the dynamics occurs as for the Lang-Kobayashi model at a large number of delays.

FIG. 11 .
FIG. 11.(a) Correlation time and (b) permutation entropy for Mackey-Glass versus the number of delays.As in the attractor plots for M < 705, high complexity can be observed.A sudden complexity collapse occurs at M = 705.

TABLE I .
Parameter values used for the simulation of the Lang-Kobayashi model Eq.(1).