Dissipative systems with nonlocal delayed feedback control

We present a linear model, which mimics the response of a spatially extended dissipative medium to a distant perturbation, and investigate its dynamics under delayed feedback control. The time a perturbation needs to propagate to a measurement point is captured by an inherent delay time (or latency). A detailed linear stability analysis demonstrates that a non-zero system delay acts destabilizing on the otherwise stable fixed point for sufficiently large feedback strengths. The imaginary part of the dominant eigenvalue is bounded by twice the feedback strength. In the relevant parameter space it changes discontinuously along specific lines when switching between branches of eigenvalues. When the feedback control force is bounded by a sigmoid function, a supercritical Hopf bifurcation occurs at the stability-instability transition. The frequency and amplitude of the resulting limit cycles respond to parameter changes like the dominant eigenvalue. In particular, they show similar discontinuities along specific lines. These results are largely independent of the exact shape of the sigmoid function. Our findings match well with previously reported results on a feedback-induced instability of vortex diffusion in a rotationally driven Newtonian fluid [M. Zeitz, P. Gurevich, and H. Stark, Eur. Phys. J. E 38, 22 (2015)]. Thus, our model captures the essential features of nonlocal delayed feedback control in spatially extended dissipative systems.


Introduction
Many phenomena in soft matter science require exciting a dissipative material. From mixing two liquids [1], sorting colloids [2,3], controlling reaction rates [4] and heat transport in microfluidic devices [5], to fluid optics [1] and spiral patterns in liquid crystals [6], soft matter systems often display their most useful or interesting properties under external stresses. Several control and driving schemes have already been applied to these systems, including optimal control [7], hysteresis control [8], and time-delayed feedback [9,8,10]. These methods sense the characteristic response of a material and adapt their control to it. Here, we focus on time-delayed feedback control of a linear model system with internal delay, which captures the essential features of how soft matter systems respond to such a feedback scheme.
Time-delayed feedback is a control strategy that was originally proposed by Pyragas to control chaotic systems and stabilize unstable periodic orbits [11,12]. It has since been applied to various dynamical systems, such as lasers [13,14,15] and neural networks [16]. The method falls into the broader category of closed-loop control schemes because its control force depends purely on the present and past states of the controlled system. Often, time-delayed feedback is called a noninvasive control scheme, as its stabilizing force ideally vanishes in a stabilized state [17].
Earlier theoretical studies [18,8,10] and experiments [9,19] applied delayed feedback to the spatiotemporal dynamics of specific systems. While several investigations of purely temporal systems with an intrinsic latency have focused on stabilizing unstable foci [20,21,22], in one earlier study on vortex diffusion at low Reynolds numbers an oscillating fluid flow in a circular geometry was initiated by delayed feedback [8]. Thus, time-delayed feedback used invasively can also destabilize stable fixed points and thereby potentially create novel nonequilibrium states in soft matter systems.
In spatially extended systems the response to a control force at location F needs the system or intrinsic delay time to propagate to a distant location A. There, it is sensed and then fed back into the system at the first location F (see the schematic in Fig. 1). To mimic this situation, we propose a linear dissipative model system, with the inherent system delay time appearing in the response function. We then apply additional time-delayed feedback control and concentrate on the case, where the system remains stable for vanishing system delay. We perform a detailed linear stability analysis and demonstrate how a nonzero system delay acts destabilizing on the otherwise stable fixed point. Typically, the response of a physical system to an external stimulus is bounded by nonlinearities. To mimic such a behavior, we introduce a sigmoid function to limit the strength of the feedback control force. As a result, stable limit cycles appear in the unstable parameter regions. While their amplitudes depend on the specific realization of the sigmoid as hyperbolic tangent, algebraic sigmoid, and ramp function, their frequencies are similar. Both, the stability-instability transition and the appearance of limit cycles, match to the findings in Ref. [8]. This suggests that our model captures the essential features of a spatially extended dissipative systems when subjected to nonlocal delayed feedback.
We present the linear dissipative model in Sec. 2, derive the characteristic function of its fixed point, and describe our numerical methods. In Sec. 3 we investigate and discuss the linear stability of the fixed point. In Sec. 4 we modify the feedback term by the sigmoid function such that its absolute value is limited and study the stable limit cycles arising from this modification. Finally, we summarize our findings and conclude in Sec. 5. Figure 1: Schematic of the model with nonlocal delayed feedback. The process A(t) measured at location A feeds back to the time-delayed control force F (t) at location, the response of which then propagates during time τ s to location A.

Derivation
We derive a simple model for a dissipative physical process A(t), which is driven by an external, nonlocal force F (t). Generally, the linear response of A(t) to F (t) is written as where the response function χ(t) characterizes the physical system completely. If F acts at some distance from the position where A is observed, there will be a system delay (or latency) quantified by time τ s , before A is affected by F . We describe this causal link in χ using the Heaviside function Θ(t − τ s ). Furthermore, in dissipative systems the response to some perturbation often decays exponentially with rate α and the response function becomes Substituting χ(t) into Eq. (1) and differentiating both sides with respect to t, we find Following Pyragas [11], we now implement for F (t) delayed feedback control with delay time τ c , control strength κ, and a constant external force F 0 , and substitute the expression into Eq. (3): This equation has one fixed point A * , where ∂A/∂t = 0, which is determined by the first two terms because the delayed control term vanishes for constant solutions, We nondimensionalize time t and κ with α, force F 0 with αA * , and A with A * to obtain In the following we study this form of the delay differential equation (DDE). Note that for κ < 0 the delayed feedback term acts to destabilize A, because when the feedback term gives a positive contribution to ∂A/∂t. Already for τ s = 0, this destabilizes the system for sufficiently large |κ|. In the following we take κ ≥ 0 to explore the role of τ s in the model.

Characteristic function
Within control theory Eq. (7) is categorized as a closed loop system with delayed feedback. Because such systems can become unstable [23], we investigate the linear stability of the fixed point at A = 1 by introducing the perturbation ansatz with |ε| 1 and complex eigenvalues λ ∈ C into Eq. (7), which gives a transcendental equation for λ, In the following we solve this equation numerically and study the properties of its solutions. The eigenvalues λ are roots of the characteristic function g(λ), Note that they cannot be expressed using the Lambert W function, as is usually possible for time-delay differential equations [17], due to the double delay terms. However, we can infer some properties of the roots from the structure of g. First, with each complex λ also its complex conjugateλ is a root of g(λ) since κ, τ s , and τ c are real parameters. Second, in appendix B we show that any root λ of g with Re(λ) > 0 has a nonzero imaginary part bounded by 2κ, |Im(λ)| < 2κ. Thus any eigenvalue associated with an unstable fixed point has a nonzero imaginary part, which is bounded by control strength. Since the imaginary part is the oscillation frequency, we conclude that fast oscillations require sufficiently large control strengths.

Root finding
In general, the roots of the characteristic function g(λ) cannot be determined analytically. Therefore, we use a numerical root finding algorithm to find the dominant eigenvalue in our stability analysis, i.e., the eigenvalue with the largest real part for one set of system parameters. As a prerequisite we need to restrict the complex plane to a finite box, which is guaranteed to contain the dominant eigenvalue. In appendix B we prove the dominant eigenvalue for any parameter combination is contained in the compact complex region B ⊂ C.
We search a rectangular superset of B (see appendix) using an interval Newton method, which provably finds all complex roots of g(λ) in B. For example, it is described in Ref. [24] and implemented in Ref. [25].

Trajectories
We also calculate the time evolution of A(t) to investigate its long-time behavior, such as limit cycles in the case of a bounded control force. To do so, we use the method of steps based on a 5th-order Runge-Kutta method [26] and the 4th-order Rosenbrock method RODAS [27] when using the ramp function to bound the control force (see Sec. 4). Both methods are implemented in the software package DifferentialEquations.jl [28,29,30]. The method of steps treats delay terms by splitting the domain of integration into a sequence of time intervals, so that on each interval the delayed values of the dynamic quantities are fully known. In our case [see Eq. (15) in Sec. 4] the method initially requires a history function describing are then known until t = τ s . For remaining intervals nτ s < t ≤ (n + 1)τ s (with n ∈ N) integration continues using A(t) from the previous intervals.
All numerical calculations are performed using the Julia programming language [31,32] and all plots are created using the matplotlib package [33].

Stability for vanishing delays
In the limiting case of vanishing system delay, τ s → 0, the fixed point is always stable. To prove this, we set τ s = 0 and take the real part of g: We prove by contradiction: Suppose that Re(λ) ≥ 0 together with κ, τ c > 0. Then the final term in Eq. (12) is the only (potentially) negative contribution necessary to obtain Re(g) = 0. However, because the absolute value of this term is always smaller than or equal to κ, roots with Re(λ) ≥ 0 cannot exist. Therefore, all eigenvalues have Re(λ) < 0 and the fixed point must be stable for τ s → 0.
In the limit of vanishing control delay, τ c → 0, the fixed point is stable because the characteristic function lim τc→0 g(λ) = λ + 1 (13) only has one (real) root at λ = −1. This is also clear from Eq. (7) since the delay term vanishes completely.

Stability-to-instability transition
A fixed point is unstable if Re(λ) > 0 for its dominant eigenvalue λ. Based on numerical calculations of the dominant eigenvalues over a wide set of parameter combinations, we make several observations, which we report here. For each τ s > 0 and τ c > 0, there is a critical control strength κ crit (τ s , τ c ) determined by Re(λ) = 0 beyond which the fixed point is unstable for all κ > κ crit . The κ crit form a manifold in parameter space; cross sections for various τ s are displayed in Fig. 2. For τ c → 0 the critical value κ crit diverges like τ c −1 . In this limit g(λ) is approximated to linear order in τ c as from which follows that any root λ stays the same as long as the product κτ c remains constant. For increasing τ c specific values of κ crit exist, where the unstable eigenvalue switches from one branch to another. We observe the resulting cusps become less prominent as τ c increases and eventually κ crit as a function of τ c approaches a constant value for each τ s . Notably, these cusps imply that for κ close to the cusp value alternating regions of stability and instability occur as τ c increases. We also observe that for small Since also Re(λ) = 0 at the transition, there is a supercritical Hopf bifurcation. Furthermore, as τ c increases the imaginary part repeatedly decreases and then jumps to a larger value when the dominant eigenvalue switches to another branch. Figure 3(a) shows that the jump lines correspond to valleys in the real part of the dominant eigenvalue. We furthermore note the similarities between these dominant eigenvalues and eigenvalues described in Ref. [8] for the specific case of delayed feedback control applied to vortex diffusion in a circular geometry. The characteristic function for that system is derived by solving the spatiotemporal problem explicitly. It contains Bessel functions, which play a similar role as the exponentials in Eq. (10). Furthermore, the diffusive response function [compare Eq. (17) in appendix A with α → 0], has an initial increase and then decays to zero, just as in our model. We consider these similarities a strong indication that some spatiotemporal systems map well on our simplified model response functions approximately given by Eq. (2).

Bounded control force
Typically, the response of a physical system to an external stimulus is bounded by nonlinearities. We mimic this behavior here by implementing a bounded control force. Bounded delayed feedback was previously studied, e.g. to suppress overshooting due to overcompensating control forces [34], or for hydrodynamic vortex diffusion in a circular geometry to stabilize unstable modes [8].

Model and setup
We start by modifying the DDE of Eq. (7) such that the time-delayed control term is limited by a monotonically increasing odd function σ(x) with σ(±∞) = ±1 and σ (0) = 1: Compared to Eq. (7) we set the constant force to zero, which shifts the fixed point of the modified DDE with the bounded feedback to A * = 0. Linearizing around this fixed point, one recovers Eq. (7) without the term +1.
We solve Eq. (15) numerically as described in Sec.

with the history function
A(t < 0) = 0 and a small initial perturbation, A(t = 0) = 10 −3 . 1 If the fixed point is unstable, the perturbation will grow over time, otherwise it will decay to zero. Using this setup, we study the long-time behavior of the bounded system and its relationship to the unstable fixed point studied in Sec. 2 by calculating the time evolution of A(t) until time T = 10 3 and by examining the frequencies and amplitudes of the occurring limit cycles for times 0.8 T < t < T .

Long-time dynamics
For all three sigmoid functions σ, the bounded system displays a Hopf bifurcation, where the fixed point A(t) = 0 becomes unstable and a limit cycle emerges (see Fig. 5). It is stabilized by the upper and lower bound of the control force, which would otherwise grow to infinity. The limit cycle does not correspond to an unstable periodic orbit of the uncontrolled system as envisaged in Pyragas' original idea [11], because the feedback term does not vanish when the limit cycle is reached. Generally, the limit cycles for all sigmoid functions should converge for κ → ∞ because in this limit they all approach the discontinuous step function. In the studied parameter region, the frequencies of the observed limit cycles are determined by the imaginary part of the dominant eigenvalue at the unstable fixed point. This becomes evident when comparing the frequency plots in the left column of Fig. 5 with Fig. 3. In particular, the jumps from one branch to the other agree well for both frequencies. However, the amplitude ∆A of the limit cycle  Fig. 3(a)]. It rather behaves like the frequency of the limit cycle with one difference. It increases as τ c grows and drops to a smaller value when the limit cycle frequency jumps to another branch. So, at the discontinuity lines the sudden increase in frequency is accompanied by a sudden decrease in amplitude. This is explicitly demonstrated in  There are some notable differences between the limit cycles generated by the three sigmoid functions bounding the control term. Most obviously, their amplitudes (right column of Fig. 5) behave differently at the bifurcation line: while they increase smoothly from zero for the hyperbolic tangent and algebraic sigmoid functions as in a supercritical Hopf bifurcation [rows (b) and (c)], they jump to a nonzero value for the ramp sigmoid function [row (a)]. The ramp function is a special case because it is linear up to the bounding values. This causes the control amplitude to always reach its maximum value in the unstable region, once the transition line is crossed. The step-like behavior could, for example, be used in experiments to accurately locate the transition line. Furthermore, for the ramp function the amplitudes of the limit cycles are largest close to the bifurcation line, while they increase for the smooth functions when moving away from the bifurcation line with increasing κ. Finally, a closer inspection shows that the discontinuity lines of the limit-cycle frequencies for the ramp function accurately track the corresponding lines in the imaginary part of the dominant eigenvalue in Fig. 3(b). However, there are slight deviations for the smooth sigmoid functions.

Transient pulse trains
In particular, for large control times τ c we observe transient pulse trains at the beginning of the dynamic evolution of our system. They occur for stable and unstable fixed points with decaying or growing amplitude, respectively. One example, when the fixed point is unstable, is displayed in Fig. 7(a). These pulse trains grow into regular limit cycles, as shown in Fig. 7(b), while for stable fixed points their amplitude decays to zero. Generally, we observe that pulse trains repeat with a periodicity given by τ c and their oscillation frequency is close to the imaginary part of the dominant eigenvalue Im(λ)/2π.

Conclusions
Motivated by the growing interest in applying delayed feedback to soft matter systems, we have studied a linear dissipative model, which mimics nonlocal delayed feedback coupling. To do so, we introduced an inherent system delay in the response function, which represents the time a perturbation needs to propagate to a distant location. Our results provide an indication how nonlocal delayed feedback in a spatially extended system determines its dynamics and helps to classify the observed spatiotemporal response.
In our investigations we concentrated on the case where the system remains in its stable fixed point when the intrinsic delay is not present. Turning on the intrinsic delay, the feedback-driven system becomes unstable for sufficiently large control strengths. We are able to show that the absolute imaginary part of the dominant eigenvalue is bounded from above by twice the feedback strength. Stability-instability transition curves and the imaginary part of the dominant eigenvalue match well with the findings reported in Ref. [8], where delayed feedback was applied to vorticity diffusion of a Newtonian fluid in a circular geometry at low Reynolds numbers. This demonstrates that our simple model system captures the essential properties of the spatiotemporal dynamics in a specific system.
To further study the feedback-induced instability, we examined the long-time dynamics of our model with bounded feedback, which in linearized form yields the original system. Realizing the bounded feedback by smooth sigmoid functions, the stabilityinstability transition occurs via a supercritical Hopf bifurcation. Thus, the fixed point becomes unstable and a stable limit cycle evolves continuously. In contrast to Pyragas' original idea [11], this limit cycle does not correspond to an unstable periodic orbit of the uncontrolled system, but is stabilized by the control force. For all three sigmoid functions and across many parameter combinations, the frequencies of the limit cycle match well the imaginary part of the fixed point's dominant eigenvalue. Discontinuity lines are visible, which occur when the dominant eigenvalue switches from one branch to another. When crossing these lines by a small increase in control delay, the limit-cycle frequency jumps up and the amplitude drops sharply. The results for the nonsmooth ramp function differ from the other sigmoid functions, because at the Hopf bifurcation the amplitude of the limit cycle jumps to a nonzero value. In addition, the discontinuity lines of the limit-cycle frequencies accurately track the corresponding lines for the imaginary part of the dominant eigenvalue of the linearized system. Both features predestine the ramp sigmoid function to accurately determine the stability-instability transition in an experimental system. The linear dissipative model presented in this article with its intrinsic delay time helps to classify and understand the spatiotemporal response of a spatially extended system subject to nonlocal delayed feedback. Our work demonstrates that this model already exhibits complex and non-trivial behavior. It also provides an example for double delay systems, which have recently attracted attention [36,37,38]. In future studies we will apply delayed feedback to specific nonlinear soft matter systems such as photoresponsive fluid interfaces [39] and viscoelastic flow in Taylor-Couette geometry [40]. The work presented in this article will help us to categorize the observed spatiotemporal dynamics.
where D is the diffusion constant. As initial condition at t = 0 we choose a delta peak located at r = 0, ρ(r, 0) = N δ(r). The solution of this problem is given by where r = |r| is the distance from the perturbation. The density at any point r 0 = 0 increases up to the time and then decreases to zero. Note that for pure diffusion (α → 0) the maximum is reached at t 0 = r 2 0 (4D) −1 . Thus, a disturbance initiated at r = 0 needs the intrinsic delay time t 0 to reach r 0 . To approximate this behavior in a response function of the form given in Eq. (2), we assume a step function where the substance ρ jumps to its maximum value ρ(r 0 , t 0 ) and then decays exponentially, Thus, τ s = t 0 is the intrinsic delay time and the physical decay rate α of our substance enters directly the response function.

B Search region for dominant eigenvalues
We find the roots of the characteristic function g numerically. Because a numerical search on the infinite complex plane is unfeasible, we restrict our search to a bounded region which is guaranteed to contain the dominant eigenvalue, i.e., the root of g with the largest real part. First, we find an upper bound to the real part of all eigenvalues using Re[g(λ)] = 0: The cosines may at most change the signs of the terms such that both contribute positively. Thus, Assuming Re(λ) ≥ 0, and we solve for Re(λ) using Lambert's W function, With our previous assumption we have Notably the upper bound is independent of τ c . Second, we find a lower bound for Re(λ) such that our search region contains at least one eigenvalue. To simplify our search, we concentrate on λ ∈ R, which is determined by On the interval (−∞, 0] the l.h.s. of this equation goes continuously from −∞ to 0 and the r.h.s. goes continuously from +∞ to −1. Because both sides are continuous and strictly monotonic for τ c , τ s > 0, they share exactly one value in the overlap, i.e. −1 ≤ λ ≤ 0. Therefore, a search region with Re(λ) ≥ −1, which also contains the real axis, will always contain at least one (real) eigenvalue. Third, we find an upper bound for the imaginary parts of all eigenvalues using Im[g(λ)] = 0. Im(λ) = κ exp[−τ s Re(λ)] sin[τ s Im(λ)] − κ exp[−(τ s + τ c )Re(λ)] sin[(τ s + τ c )Im(λ)] (26) Here, we simply drop the sines, assuming all terms contribute positively, The r.h.s. depends on Re(λ) and on the interval determined for Re(λ). It is maximal for Re(λ) = −1 implying Im(λ) ≤ κ e τs [1 + e τc ]. Note that the symmetry g(λ) =ḡ(λ) implies that we only need to search the positive half of the complex plane, i.e., Im(λ) ≥ 0, because all complex roots appear in conjugate pairs. In summary, we have shown that the following set of complex numbers B must always contain the eigenvalue with the largest real part: Both regions are displayed in Fig. 8. Note that the shape of B immediately implies that any eigenvalue λ with Re(λ) > 0 has |Im(λ)| < 2κ. Furthermore, it can be shown that Re(λ) > 0 implies |Im(λ)| = 0 because Eq. (25) has no solutions for λ > 0. In this case the l.h.s. of Eq. (25) is always positive while the r.h.s. is always smaller than −1 and therefore no solution with Re(λ) > 0 and |Im(λ)| = 0 exists.