Delayed feedback control of unstable steady states with high-frequency modulation of the delay

We analyze the stabilization of unstable steady states by delayed feedback control with a periodic time-varying delay in the regime of a high-frequency modulation of the delay. The average effect of the delayed feedback term in the control force is equivalent to a distributed delay in the interval of the modulation, and the obtained distribution depends on the type of the modulation. In our analysis we use a simple generic normal form of an unstable focus, and investigate the effects of phase-dependent coupling and the influence of the control loop latency on the controllability. In addition, we have explored the influence of the modulation of the delays in multiple delay feedback schemes consisting of two independent delay lines of Pyragas type. A main advantage of the variable delay is the considerably larger domain of stabilization in parameter space.


I. INTRODUCTION
The possibility to stabilize unstable periodic or stationary states embedded in chaotic attractors has been elaborated more than two decades ago [1]. The main idea consists of using small external perturbations to force the system to follow one of its stable manifolds. These perturbations are applied at specific instances when the chaotic trajectory is close to the desired periodic orbit. This seemingly straightforward theoretical concept has caused a revolution in applied nonlinear science. It has been realized that the control of chaos could have a significant outcome in real-world experiments, where one could generate different kinds of ordered behavior from an utterly erratic one [2].
A conceptually simple method to stabilize unstable equilibria and periodic orbits is the time-delayed feedback control (TDFC) introduced by Pyragas [3,4]. Here, the perturbation has the form of a continuous feedback constructed from the difference between some suitable scalar signal obtained from the system and the same signal delayed by a constant time τ . The difference signal is amplified by the gain factor K and then re-injected into the original system. For a certain choice of the feedback gain K and the delay time τ , the control of the unstable state can be realized, in which case the feedback force vanishes by construction, making the method essentially noninvasive.
Since a detailed knowledge of the target state is not required and the controller is very robust with respect to * Electronic address: agjurcin@pmf.ukim.mk † Electronic address: schoell@physik.tu-berlin.de noise, the Pyragas control has become one of the most popular control methods used in experiments and even technological applications. The application is quite diverse, and includes stabilization of unstable states in electronic chaotic oscillators [5][6][7], mechanical pendulums [8,9], laser systems [10][11][12], electrochemical systems [13,14], drift waves in a magnetized laboratory plasma [15], chaotic Taylor-Couette flow [16], cardiac systems [17,18], ferromagnetic resonance systems [19], gas discharge systems [20,21], controlling helicopter rotor blades [22], controlling the walking mechanism of a robot [23][24][25], stabilization of cantilever oscillations in an atomic force microscope [26], and controlling librational motion of a tethered satellite system in an elliptic orbit [27], amongst others. For a comprehensive review of the technical implementation of the method, see [2,4]. In parallel to experimental realization, substantial work has been done to understand the control mechanism analytically [28][29][30][31]. A notable result has been reported recently in the context of refuting the alleged oddnumber limitation [32,33], believed to be a severe limitation of the delayed feedback control technique for almost a decade [34]. A way to overcome the odd-number limitation was proposed by Pyragas by using an additional unstable degree of freedom in a feedback loop [35][36][37].
The original Pyragas method was subsequently improved by introducing multiple delayed signals into the feedback control force, both with commensurate or incommensurate delay times [38][39][40][41]. Another improvement of the method was achieved by introducing a timevarying delay into the delayed argument of the control force which could be realized experimentally by changing some characteristics of the delay line in a deterministic (periodic) or a stochastic fashion [42][43][44]. For example, in coupled laser systems the delay time could be modu-lated by periodically changing the distance of the lasers or external resonator, e.g. by piezoelectric modulation. For systems acting on a slower timescale, e.g. electronic or mechanic devices, one may take advantage of digital delay lines instead, where the time-functional form of the delay modulation could be controlled by an external clock frequency modulator [45]. In this way, one could practically realize the variable-delay feedback controller and stabilize unstable steady states and periodic orbits over a much larger domain of control parameters. To keep the method noninvasive also in case of unstable periodic orbits, the delay modulation should be chosen appropriately, e.g., in a form of a (non)periodic change of the delay between multiples of the basic period of the uncontrolled unstable orbit [46].
In this paper, we aim to provide a deeper analytical insight into the mechanism of stabilization of unstable equilibria by the delayed feedback method with a timedependent delay time, and extend this concept by considering the influence of a phase-dependent coupling, the control loop latency, and multiple delay lines on the controllability, which is relevant in real experiments. A general description of the variable-delay feedback controller is given in Sec. II, where it is shown that the method could be explored analytically in the domain of a highfrequency modulation of the delay. In this case, we use the formalism of distributed delays [47,48], where the average contribution of the time-varying delay is represented by an integral kernel describing a particular delay distribution. In Sec. III, we consider a phase-dependent coupling of the control force, which is a particular nondiagonal generalization of the standard diagonal coupling scheme. As a model subject to control we use a generic two-dimensional linear system that has an unstable steady state of focus type. This model mimics the two-dimensional center manifold of a general nonlinear system, capturing the essence of the dynamics in the vicinity of its unstable fixed point. In Sec. IV, we investigate the influence of control-loop latency on the efficiency of the controller. In Sec. V, we extend our method to the case of multiple delay feedback by considering two independent delay lines with time-varying delays. We conclude our findings in Sec. VI.

II. VARIABLE-DELAY FEEDBACK CONTROL
We consider a general n-dimensional dynamical systeṁ x = F(x), where F is the vector field describing the system's dynamics, and x = x(t) ∈ R n is the state column vector. The equilibria x * i of the system are the solutions of F(x * i ) = 0, and the stability of a particular equilibrium x * is determined by the eigenvalues of the Jacobian matrix A = ∂ x F(x) calculated at x = x * . In the following we assume that x * is an unstable equilibrium, having at least one eigenvalue with a positive real part.
Under a linear variable-delay feedback control, the original system is transformed intȯ where K is an n × n feedback gain matrix, and τ (t) is the time-dependent delay time. By choosing K and τ (t) appropriately, we aim to stabilize the unstable equilibrium x * . Without loss of generality, we assume that the equilibrium has already been moved to the origin by a translation of axes, such that x * = 0. We thus linearize the controlled system in the neighborhood of the origin to obtainẋ We consider a deterministic modulation of the variable time-delay τ (t) in a fixed interval around a nominal (average) delay value τ 0 where f : R → [−1, 1] is a 2π-periodic function with zero mean, and ε and ν are additional variable-delay control parameters denoting the amplitude and the (angular) frequency of the modulation, respectively. The stability of the origin can be inferred by numerically integrating the linear variable-delay system (2) for different values of the control parameters K, τ 0 , ε and ν, thus determining the domains in the parameter space for which the stabilization can be realized.
In the limiting case of a high-frequency modulation [43,49], when the parameter ν becomes sufficiently large compared to the uncontrolled system's dynamics described by an intrinsic frequency ω, the linearized system (2) with a modulated delay has the same asymptotic stability properties as the averaged linear distributed-delay systeṁ The probability density function ρ(θ), i.e. the distributed-delay kernel, is defined in a way that ρ(θ)dθ gives the fraction of time for which τ (t) lies between θ and θ + dθ, satisfying ρ(θ) ≥ 0 and the probability normalization condition In practice, the transition to the distributed-delay limit case does not require very large modulation frequencies, and therefore variations of the delay in the experiment are a very practical way to create different types of delay distributions. For a continuous modulation in an ε interval around τ 0 in form of a sawtooth-wave (triangular) modulation of the delay, ρ(θ) is a constant in the interval of the delay variation. Under the probability normalization condition (5) we obtain a uniform distribution which does not depend on the skewness of the sawtooth function. For a sine-wave modulation, the time interval dt during which the delay τ changes by dτ is given by The fraction of time dt within a half-period π/ν of the delay function τ is equivalent to the product ρ(τ )dτ . In terms of our previous notation, we get A periodic change of τ (t) between two fixed values τ 0 − ε and τ 0 + ε for the same time duration, results in a square wave (rectangular) modulation with a two-peak distribution In an analogous way, we may obtain the corresponding distributed delay kernels for other types of delay modulations.
Applying the exponential ansatz x(t) ∼ exp(Λt) in Eq. (4), we obtain the characteristic equation for the eigenvalues Λ (13) where I is the identity matrix. Since ρ(θ) is nonzero only between τ 0 − ε and τ 0 + ε, we have The quantity summarizes the effect of a given modulation, and for the above modulation types reads , sawtooth-wave, I 0 (Λε), sine-wave, cosh(Λε), square-wave, in the vicinity of a steady state (red/gray horizontal line) of a system subjected to variable-delay feedback control with a fast sawtooth wave. The corresponding uniform delay distribution (7) creates a sliding average between t − τ0 − ε and t − τ0 + ε (area highlighted in blue/gray) which at time t results in the reference signal (blue/gray, oscillating). The control force is constructed from the difference signal as marked by the vertical black arrow.
where I 0 is the modified Bessel function of the first kind of order 0. In the non-modulated case, χ(Λ, 0) ≡ 1, we have the usual TDFC. For general distributed delay, the delay term exp(−Λτ 0 ) is replaced in the characteristic equation (13) by the Laplace transform of the delay kernel ρ(θ) [47].
The function χ(Λ, ε) also allows for a compact explanation of the mechanism of variable-delay feedback control (VDFC). The expressions listed in Eq. (16) have regions {Λε ∈ C} for which |χ| < 1, mostly for Re(Λ) ≪ Im(Λ). In particular, values of χ ≈ 0 can occur, meaning that the delay term effectively has been suppressed and vanishes from the characteristic Eq. (14). In the original equations of motion (4) this situation can be interpreted as destructive interference of the delay signals covered by the distribution. Positive and negative phases of the spiral oscillation in the neighborhood of the steady state cancel each other out creating a reference term x ref (t) ≈ x * for a simplified controlleṙ The same idea can be regarded as the motivation for the original TDFC method. The advantage of VDFC lies in an improved approximation of the target state by the delay terms. Fig. 1 illustrates the control mechanism for an almost ideal situation, in which the steady state is approximated by the averaged delay signal. The instantaneous part of the coupling term stabilizes the system while pulling it towards the reference signal x ref (t) which due to the small χ-value has a smaller amplitude than the remaining instantaneous x(t). The resulting stabilization of the steady state shows a high robustness against parameter detuning, because the simplified con-trol mechanism is insensitive to the phase relation between actual signal and reference signal and also works for a wide range in the coupling gain K. However, the described scenario relies on a small χ-value, which in general is not trivial to obtain. The delay terms can almost never be suppressed completely, so that only a rigorous consideration of the characteristic equation (14) reveals the full capability of the control method. Eq. (14) is transcendental with respect to Λ, possessing an infinite set of complex solutions {Λ i } ∈ C. The origin can be stabilized for those values of the control parameters ( K, τ 0 , ε) for which all the eigenvalues {Λ i } have negative real parts. The stability domain in the parameter space ( K, τ 0 , ε) can be calculated numerically given the modulation distribution. Note that the control parameter ν is lost in the transition to distributed delays. As in the usual TDFC and extended TDFC control schemes, the presence of torsion is necessary for the proposed control method to be able to stabilize equilibria. Torsion means that in its unstable subspace the steady state is only of the spiral type. To show this property, we consider the characteristic quasipolynomial H(Λ) = det Λ I − A + 1 − e −Λτ0 χ(Λ, ε) K . It is easily shown that this quasipolynomial is positive for Λ → ∞, and for , where e i are the eigenvalues of A. If A possesses an odd number of positive real eigenvalues, then H(0) < 0, and there exists at least one positive real root of H(Λ) = 0, meaning that the fixed point cannot be stabilized by the proposed method. Although the analysis is done in the distributeddelay limit, numerical simulations show that this oddnumber limitation persists also for low-frequency modulations of the delay τ (t).
The results of this section will be exploited in the following to perform a comparative analysis between the variable-delay feedback control in a distributed-delay limit and the standard delayed feedback control by using a simple normal form model as a representative of a large class of nonlinear dynamical systems. Specifically, we will investigate the effects of including a phasedependent coupling matrix and the influence of nonzero latency times in the control scheme.

III. PHASE-DEPENDENT COUPLING
We will apply the variable-delay feedback control to an unstable steady state of focus type. This system represents a generic model of an unstable steady state slightly above a Hopf bifurcation. In center manifold coordinates, the dynamics of the system is given by Eq. (2), where x = (x, y) T is the two-dimensional column vector, and A is a 2 × 2 matrix that describes the dynamics in the absence of control.
Since the stability of the free-running system is determined by the eigenvalues λ ± i ω of A, choosing λ > 0 and ω = 0 we model an unstable focus at the origin.
In this section, we investigate a phase-dependent coupling of the control force, which is relevant, e.g., in stabilization of laser devices where the optical phase occurs as an additional degree of freedom [50,51]. Such couplings have also been used in the context of refuting the odd-number limitation of delayed feedback control in autonomous systems [32,33,52], to anticipate chaos synchronization [53], and more recently, to control different synchronous oscillatory states of delay-coupled oscillator networks [54]. The phase-dependent coupling is realized via a rotational matrix K containing a variable phase ϕ, and it enters as an additional multiplicative factor to the control force where K is the feedback gain. In the distributed-delay limit, the system becomes and the stability of the origin is determined by the roots Λ of the characteristic equation This equation can be further simplified to The control method is successful if there exists a nonempty set in the four-dimensional parameter space (K, τ 0 , ε, ϕ) for which the real parts of all the character- The grayscale (color code) depicts only those values of the control parameters for which the largest real part of the complex eigenvalues Λ is negative, indicating a successful control. The depicted solid lines that enter into a description of the stability boundary are calculated from Eqs. (26) and (27), and the dashed lines from Eqs. (26) and (28). For clarity of the picture, we depict only a few branches of the boundary curves.
istic roots Λ are negative. The domain of control in the plane parametrized by the feedback phase ϕ and the feedback gain K is given in Figs. 2 and 3. The grayscale (color code) indicates only those control parameters (ϕ, K) for which the leading eigenvalues have negative real parts, and the control is more robust as these values are more negative. Panels The chosen nominal delay τ 0 = 1 is an optimal value in the non-modulated control case and ϕ = 0 (diagonal coupling). From the corresponding control domains, it is observed that an increase of the modulation amplitude ε leads to an enlargement of the control domain in the direction of the positive ϕ-axis, thus destroying the symmetry of the stability island centered at ϕ = 0. A similar result is observed for a sine-wave modulation in panels (a)-(b) in Fig. 3, corresponding to ε = 0.5 (a) and ε = 0.9 (b). A different behavior is observed in panels (c)-(d) in Fig. 3 depicting the corresponding stability domains for a square-wave modulation. In the latter case, the stability domain increases more rapidly as ε is increased from zero, until ε reaches a critical value after which the domain decreases, resulting in an instability stripe at ϕ = 0 when ε = 1. This non-monotonic change of the stability domain is typical for a square-wave modulation, and it has been observed in other circumstances.
In order to understand the mechanisms behind these numerical results, we first recall the intuitive argument about destructive interference and the simplified controller from the previous section. From Eq. (20) we directly see that the stabilizing diagonal elements of the instantaneous control term are weighted with cos ϕ. Since the onset of stability K min is dominated by the competition between λ > 0 and −K cos ϕ < 0, we expect K min ∝ λ/ cos ϕ, see also the detailed reasoning in Eq. (30) below. For |ϕ| > π/2 the instantaneous part of the coupling becomes repelling. Therefore with regard to the depicted control mechanism we expect the control to largely fail in this region. All simulations clearly support this picture. Small exceptional regions of stability exceeding this boundary arise from the non-vanishing delay terms which can support stabilization given an optimal parameter constellation. However, for large values of the mean delay τ 0 this type of support fails and |ϕ| = π/2 becomes a strict limit of control.
The asymmetry can be explained by the interaction of the non-diagonal elements of system terms and control terms in Eq. (20). If ω were zero, the control term would have no preference for positive or negative values of ϕ because the change of ϕ → −ϕ would be equivalent to mirroring the complete system at one coordinate axis, e.g. x → −x. A non-zero value of the spiral frequency ω breaks the symmetry and leads to different resulting frequencies Ω = Im(Λ) under control for different signs of ϕ. From the interference mechanism one can conclude, that at least for the continuous distributions (7) and (10) high frequencies are more favorable for successful control, because a low χ-value corresponds to many oscillation periods covered by the distribution. Indeed, if we have a closer look at the parameter constellation underlying Fig. 2, we see that for negative ϕ the imaginary part of the most unstable mode tends to vanish, thus inhibiting the interference effect and leading to control failure. In contrast, for positive ϕ the frequency increases with increasing K, thus improving control efficiency. A parametric plot of different Λ i in the complex plane reveals this property, see Fig. 4.
A detailed analytical investigation of the control domain can be done by using the characteristic equation (22) and noticing that the transition from instability to stability occurs at the boundaries of the control domain, where the leading eigenvalues Λ are purely imaginary.
Setting Λ = iΩ in Eq. (22), and separating real and imaginary parts results in a system of two real-valued equations (24) which is an implicit parametric representation of the boundaries of the control domain in the parameter space (K, τ 0 , ε, ϕ), in which the eigenfrequency Ω has the role of a parametrization variable. We have taken into account that for the above mentioned distribution kernels, ρ(θ) is even around the nominal delay τ 0 , meaning that χ(iΩ, ε) is a real-valued function sine-wave, cos (Ωε), square-wave, (25) where J 0 is the Bessel function of the first kind of order 0. Considering the positivity of the nominal delay τ 0 and the feedback gain K, as well as the positivity of the parameters ϕ, λ and ω used in the current analysis, elimination of the phase parameter ϕ from the last system of transcendental equations leads to an expression for K in terms of the eigenfrequency Ω Taking into account the multivalued properties of the arcsine function, one obtains in a similar manner the analytical expressions for the phase parameter ϕ in dependence on Ω where n is a non-negative integer. Equations (26)- (28) describe the boundary of the control domain in Figs. 2 and 3, where the solid and dashed lines represent branches of ϕ 1 (Ω) and ϕ 2 (Ω), respectively.
In the case of zero feedback phase (ϕ = 0), the phasedependent feedback force is reduced to a diagonal coupling. In the absence of modulation (TDFC), the domain of control in the plane spanned by the feedback gain and the nominal delay consists of several distinct stability islands centered at odd values of τ 0 , separated by regions corresponding to even τ 0 for which the stabilization fails for any K. This observation can be inferred from Eqs. (23)-(24) setting ϕ = 0 and ε = 0. When Ωτ 0 = (2n+1)π, then K min = λ/2 and Ω = ω. On the other hand, when Ωτ 0 = 2nπ, the control fails for any finite value of K. A detailed analysis of the control boundaries in this case (TDFC) is provided in Ref. [29]. By the same reasoning, in the presence of modulation (ε > 0) and zero phase, the corresponding minimum feedback gain values are , τ 0 = (2n + 1)π/ω, leading to a reconfiguration of the control domain depending on the type of the delay modulation [42]. Specifically, it is observed that the instability stripes at even values of τ 0 cease to exist, and the stability islands are starting to interconnect as soon as ε > 0. Although for a sawtooth and sine-wave modulation the enlargement of the control domain with increasing ε is monotonic, an oscillatory behavior is observed for a square-wave modulation due to the form of the χ function (i.e. cos(ωε)) in the denominator of the expressions for K min . In the latter case, the stability islands centered at τ 0 = (2n+1)π/ω are first enlarged and eventually interconnected with increasing ε, up to some ε value after which the control domain shrinks, gradually collapsing into several stability islands centered at τ 0 equal to even multiples of π/ω. For ϕ = 0 and a variable-delay control force (ε > 0), the minimum feedback gain at the above specific values of τ 0 cannot be deduced from Eqs. (23)-(24) by following the previous lines of deduction. Nevertheless, we can make approximate predictions for small values of ϕ, when Ωτ 0 is an integer multiple of π, and we have Ω ≈ ω from Eq. (24). Hence, in the regime of small values of ϕ, the minimum feedback gain is The resulting expressions for K min show that in the absence of any delay modulation and for small ϕ, the minimum feedback gain at the "optimal" delay values τ 0 = (2n + 1)π/ω is K min = λ/(2 cos ϕ), containing the dependence on ϕ via cos ϕ in the denominator. Consequently, the principal stability island in (ϕ, K) parametric plane is centered at ϕ = 0, for which K min = λ (see panel (a) in Fig. 2). On the other hand, at τ 0 = 2nπ/ω, the stabilization cannot be achieved for any K, regardless of the value of the phase parameter ϕ. For a non-zero modulation amplitude ε at small ϕ, the values of K min depend on the type and the amplitude of the delay modulation. Although the form of χ(iω, ε) for sawtooth and sinewave modulations is such that the control in the variabledelay case is generally possible for all nominal delays τ 0 ≥ ε, this is not the case for a square-wave modulation. Namely, in the latter case, χ(iω, ε) = cos(ωε), giving a non-monotonic behavior of K min with increasing modulation amplitude ε. Specifically, for ε = (2n + 1)π/ω, the control fails at τ 0 = (2n + 1)π/ω for any K, but it is optimal at τ 0 = 2nπ/ω.
To further investigate the effects of the modulation amplitude ε on the control efficiency, we have calculated the stability domains in the parametric plane (K, ε) for different modulation types and different values of the phase parameter ϕ. The results are depicted in Figs. 5, 6 and 7 for sawtooth-, sine-and square-wave modulation, respectively. Different panels in each figure correspond to different phases: (a) ϕ = 0, (b) ϕ = π/8, (c) ϕ = π/4, (d) ϕ = 3π/8. The parameters of the unperturbed system are λ = 0.1 and ω = π, and the nominal delay is fixed at τ 0 = 1 as before. As the modulation amplitude increases from zero, a considerable reconfiguration of the stability domain is observed, depending on the type of the delay modulation. Increasing the modulation amplitude enlarges the interval of K for successful control, although the enlargement is not necessarily monotonic, and the domain may consist of several disconnected intervals. The improvement of the delayed feedback controller from including time-varying delay is evident from the diagrams. For example, for a phase parameter value ϕ = 3π/8 (panel (d) in Figs. 5, 6 and 7), while the original Pyragas method fails for any K, the variable-delay feedback method is able to stabilize the fixed point for certain values of the modulation amplitude ε > 0.
In the second equation, ε is implicitly contained in χ(iΩ, ε), which is a real function in the delay modulation cases considered above (see Eq. 25). In general it is not possible to invert the function χ analytically. One should treat the last system as implicit parametric representation of the control boundaries and seek for the stability curves numerically. The boundaries of the corresponding control domains are given in Figs. 5-7 by the solid line.

IV. CONTROL-LOOP LATENCY
In an experimental realization of the control method, one has to take into account the latency of the feedback circuit due to the time required for the generation of the feedback control signal and its reinjection into the system. In laser systems, the latency is associated with the time the light takes to traverse the distance between the laser and the Fabry-Perot controller. It has been shown that latency time always acts destructively upon the control domains, reducing the effectiveness of the controller [29,30,51,55,56].
In this section, we investigate the effects of latency on the variable-delay feedback control in the distributeddelay limit. The latency time δ enters as an additional constant time-delay in the control force, and the system now reads In the distributed-delay limit, the system becomes leading to the characteristic equation that can be further simplified into  40) and (41), and the dashed lines from Eqs. (40) and (42). The parameters of the free-running system are as in Fig. 2 and ϕ = 0.
Upon separating the real and imaginary parts, one obtains a system of two real-valued equations. We again consider the control boundary given by Λ = iΩ which is given by the implicit parametric representation K cos(Ωδ) − K cos(Ω(δ + τ 0 ))χ(iΩ, ε) = λ, (37) Ω − K sin(Ωδ) + K sin(Ω(δ + τ 0 ))χ(iΩ, ε) = ±ω. (38) Figure 8 depicts the domain of control in the (τ 0 , K) parametric plane for a sawtooth-wave modulation at a fixed modulation amplitude ε = 0.5 and increasing latency time δ: (a) δ = 0.1, (b) δ = 0.2, (c) δ = 0.3, (d) δ = 0.4. The control domains for a sine-wave modulation corresponding to δ = 0.2 and δ = 0.4 are depicted in panels (a) and (b) in Fig. 9, and panels (c) and (d) are related domains for a square-wave modulation. The parameter values of the uncontrolled system are set to λ = 0.1 and ω = π as before. For a zero latency time, the control domain fills the whole depicted range of the parametric plane down to some minimum value of K, which for τ 0 equal to an integer multiple of π/ω is explicitly given by Eq. (29). It is observed that increasing latency time reduces the area and the robustness of the control domain. In analogy to Sec. II, one can derive approximate expressions for K min for small values of the latency δ , τ 0 = 2nπ/ω.
(39) In parallel to the conclusions of the previous section, inclusion of a variable time-delay in the feedback control force generally enlarges the control domain, making the control possible even for those values of τ 0 for which the control always fails at any K in the non-modulated case.
To describe the boundaries of the control domain in the (τ 0 ,K) plane analytically, we algebraically manipulate the system (37)- (38) and find two families of branches of solutions for K and τ 0 parametrized by the eigenfrequency Ω where n is a non-negative integer that takes care of the different branches due to the multivalued arcsine function involved in the boundary description. observed that increasing the modulation amplitude ε for a constant latency value generally leads to an enlargement of the K interval of successful control, depending on the type of the modulation. For a sufficiently large latency time when the Pyragas control fails (e.g, δ = 0.4 in panel (d) of Figs. 10-12), inclusion of a variable delay in the controller makes the control possible again for a suitable choice of the modulation amplitude ε > 0. Nevertheless, the control domain depends on the modulation type, and it is seen that in this case it is much less pronounced for a sawtooth-wave modulation than for the other two delay-modulation types. In a similar fashion, from Eqs. (37) and (38) one obtains the parametric representation of the stability boundaries in the (K, ε) plane in which the dependence of ε on the eigenfrequency Ω enters implicitly in the second equation via the function χ. The calculated boundaries of the control domain are the solid lines depicted in Figs. 10-12.
Overall, we see a massive destruction of the control domains with increasing latency time. Compared to the effect of a phase rotation in the previous section, the presence of control-loop latency tends to spoil the control mechanism radically. Besides the destructive interference effect, on which variable-delay feedback control mainly relies, an instantaneous feedback provides the main source of stability which normally is reflected in an extensive coverage of the parameter space with solutions of successful control. If this dissipation term is replaced by a much less effective term due to latency, the control performance is consequently lost. Thus in any experimental application of variable-delay feedback control one should split the control term if possible, so that the instantaneous part is implemented by a direct modification of the system to be controlled, and the variable delay part can be realized by separate devices.

V. MULTIPLE DELAY FEEDBACK CONTROL WITH VARIABLE TIME-DELAYS
In order to improve the control of unstable steady states, several extensions of the Pyragas method were proposed by involving multiple delay feedback terms in the control force [38][39][40][41]. A very efficient control scheme of this type was introduced by Ahlborn and Parlitz by utilizing a feedback force constructed from two or more independent Pyragas delayed feedback controllers with incommensurate delay times applied simultaneously in the control circuit [40,41]. A key result of this multipledelay extension was a successful stabilization of chaotic intensity fluctuations of a frequency doubled ND-doped yttrium aluminum garnet (ND:YAG) laser at higher pump rates, whose control was not achievable with a single delay controller.
In this section, we show that this multiple delay feedback control (MDFC) can be further improved by including time-varying delays in the associated feedback terms. As we have shown in the previous sections, the term improvement refers in the first instance to an extension of the parameter space of successful control. Under experimental conditions this feature is favorable if control parameters are drifting, cannot be adjusted precisely or the properties of the unstable equilibrium are unknown so that the optimal parameters of the constant delay control forces cannot be determined. Concerning the robustness of the control, we have observed for an optimal choice of modulation that the resulting maximum stability exponent in the stable domain is rather uniform. But since the complete response to perturbations of the controlled system is formed by the whole spectrum of exponents, the overall robustness can vary significantly more than reflected only by the maximum exponent. In particular, we have observed [45] that the excitation from the stabilized fixed point by additive noise is decreasing with increasing feedback strength in an optimal VDFC scheme although the maximum exponent does not vary much. This is another typical feature of VDFC which we regard as a major improvement compared to non-modulated methods. In the following analysis of modulated MDFC we restrict ourselves to the maximum stability exponent in order to allow for comparison with the previous sections. For the clarity of the presentation, we consider the simplest realization of MDFC consisting of two Pyragas-type delayed feedback lines with variable time delays, but the discussion can be straightforwardly generalized to more than two delay lines. The linear normal-form system reads where τ 1 (t) and τ 2 (t) are two different, time-dependent, 2π-periodic delay functions: The nominal delays are denoted by τ 01 and τ 02 , ε 1 and ε 2 are the corresponding modulation amplitudes, ν 1 and ν 2 are the modulation frequencies, and f 1 and f 2 are the modulation functions. For high-frequency modulations of the delays, the system is in the distributed delay regime: Using the exponential ansatz, we obtain the characteris-tic equation where are complex functions corresponding to different modulation types. Equation (49) can be re-cast in the simple form The parametric representation of the stability boundary is obtained by looking for solutions of Eq. (51) in the form Λ = iΩ and separating real and imaginary parts, The control parameter space is now six-dimensional (K 1 , K 2 , τ 01 , τ 02 , ε 1 , ε 2 ), but for experimental purposes it could be reduced by matching similar types of control parameters (e.g. To gain insight into the domain structure of the multiple variable-delay feedback control, we numerically analyze Eq. (51) in the parametric plane spanned by the two nominal delays τ 01 and τ 02 . In Fig. 13 we depict the resulting stability diagrams at different control parameter values obtained when the time delays in both feedback terms are modulated with sawtooth-waves of equal amplitudes (i.e. ε 1 = ε 2 at each panel). The parameters of the unstable focus are taken as λ = 0.1 and ω = π throughout this section. The stability area is given in graytones (color online), corresponding to those values of the control parameters for which the fixed point control could be achieved. Panels (a)-(d) correspond to increasing values of ε 1,2 for a symmetrical choice of the feedback gains K 1 = K 2 = 0.1. Panel (a) is the resulting stability domain without any delay variation (ε 1,2 = 0), i.e. MDFC with constant delays. The stability domain is symmetrical with respect to the diagonal τ 01 = τ 02 , and it is filled with oval-shaped instability islands (white regions) encompassing the points (τ 01 , τ 02 ) = (2nπ/ω, 2mπ/ω), where n and m are nonnegative integers. At these particular time delays, the fixed point is unstable for any values of the feedback gains K 1 and K 2 . The result could be shown analytically from Eqs. (52)-(53) by setting Ωτ 01 = 2nπ and Ωτ 02 = 2mπ. In the same manner, we derive the conditions for the minimum feedback gains at the nominal delay values between these instability points. Specifically, at (τ 01 , τ 02 ) = ((2n + 1)π/ω, 2mπ/ω) the minimum feedback gain is K min 1 = λ/2; at (τ 01 , τ 02 ) = (2nπ/ω, (2m + 1)π/ω) the minimum feedback gain is K min 2 = λ/2; at (τ 01 , τ 02 ) = ((2n + 1)π/ω, (2m + 1)π/ω) the minimum feedback gains satisfy condition K min 1 + K min 2 = λ/2. In the latter case, by tuning the value of one of the feedback gain parameters, a successful stabilization could be achieved at negative values of the second feedback gain parameter.
In panels (e)-(h) of Fig. 13 we show the corresponding increase of the control domain area for K 1 = 0.05 and K 2 = 0.2. This asymmetric choice of the feedback gains results in asymmetry of the stability domain with respect to the diagonal τ 01 = τ 02 , leading to horizontal enlargement of the instability islands in the absence of delay modulation (MDFC, panel a), which eventually connect into horizontal instability stripes at τ 02 = 2mπ/ω for lower values of the gain parameter K 1 . By including variable delays, the stability domain increases monotonically with increasing ε 1,2 , and the instability islands gradually shrink, eventually disappearing at large ε 1,2 .
A similar monotonic increase of the stability regions with increasing modulation amplitudes in both symmetric and asymmetric choice of the feedback gain parameters is also observed for sine-wave modulations of the delays, and the results are not shown for compactness of the presentation. For a square-wave modulation (Fig.  14), we observe a non-monotonic behavior of the stability area, increasing faster with ε 1,2 , attaining its maximum at ε 1,2 = 0.5 (panels a and c), and then decreasing again to the same domain structure as in the unmodulated case (ε 1,2 = 1, panels b and d), but now with the instability islands centered at the nominal delay values being odd multiples of π/ω (compare with panels a and e in Fig. 13). This oscillatory switching of the instability islands between (τ 01 , τ 02 ) = (2nπ/ω, 2mπ/ω) and (τ 01 , τ 02 ) = ((2n + 1)π/ω, (2m + 1)π/ω) continues as the modulation amplitudes ε 1,2 are further increased, and it is confirmed analytically by Eqs. (54)- (57).
It is also instructive to investigate the stability diagrams when the corresponding delays in the two feedback control terms are modulated with different modulation types. A general conclusion can be extracted from the sample of the results provided in Fig. 15. The diagrams are calculated for a sawtooth-wave variation of τ 1 (t) and a square-wave variation of τ 2 (t) at equal amplitudes (ε 1 = ε 2 ) and asymmetric gains: K 1 = 0.2 and K 2 = 0.05 in panels (a)-(b), and K 1 = 0.05 and K 2 = 0.2 in panels (c)-(d). Panels (a) and (c) correspond to ε 1,2 = 0.5 and panels (b) and (d) to ε 1,2 = 1. The overall behavior of the control domain with increasing ε 1,2 is mainly determined by the dominant modulation type, i.e., the modulation of the feedback with higher gain, and it is more pronounced as the difference between the respective feedback gain values becomes larger. By gradually increasing the modulation amplitudes, the domain is enlarged monotonically in panels (a)-(b) where the feedback term with a sawtooth-wave modulation dominates over the square-wave term (K 1 > K 2 ). The reversed case scenario (K 2 > K 1 ) in panels (c)-(d) shows a non-monotonic behavior typical for a square-wave modulation by the appearance of instability stripes oscillating between the values of τ 02 being even and odd multiples of π/ω.
To derive the equations for the boundary curves in Figs. (13)-(15), we follow the approach in Ref. [57] and rewrite the characteristic Eq. (51) as where a(Λ) and b(Λ) are given by At the control boundary (Λ = iΩ) the three terms in Eq.
(58) can be considered as three vectors in the complex plane, with the corresponding magnitudes 1, |a(iΩ)| and |b(iΩ)|. According to Eq. (58), the sum of these vectors is a zero vector, and from the triangle formed by the vectors it is straightforward to obtain the parametric representation of τ 01 and τ 02 on the eigenfrequency Ω: where u ± 0 and v ± 0 are the smallest possible integers such that the corresponding values of τ 01 and τ 02 are all nonnegative, and α 1 , α 2 ∈ [0, π] are the internal angles of the triangle formed by the vectors, calculated from the law of cosines as: To gain further insight into the superiority of MDFC with time-varying delays with respect to the fixed delays realization, we have numerically calculated the stability diagrams in the parametric plane of the nominal delay τ 01 and the feedback gain K 1 of the first feedback line at different ε 1,2 , by fixing the nominal delay τ 02 and the feedback gain K 2 of the second feedback line. For the value τ 02 we choose an even multiple of π/ω at which the fixed point control by MDFC with constant delays always fails if τ 01 is also an even multiple of π/ω. In Fig. 16 we show the corresponding control domains for τ 02 = 2π/ω = 2 and K 2 = 0.1 for sawtooth-wave modulations of the delays at different modulation amplitudes: (a) ε 1,2 =0 (MDFC), (b) ε 1,2 =0.25, (c) ε 1,2 =0.5, (d) ε 1,2 =1. For fixed delays (MDFC case, panel a), the control domain is consisted of isolated stability islands encompassing τ 01 = (2n + 1)π/ω, disconnected by instability stripes at τ 01 = 2nπ/ω at which control fails at any K 1,2 . As the modulation amplitudes increase (panels b-d), the stability islands gradually expand into a connected stability region, making the control possible at any τ 01 . This monotonic expansion of the stability domain is sustained in the case of a sine-wave modulation, but not for a square-wave modulation, as expected from the earlier analysis. In the latter case, instability stripes are disappearing and re-appearing again with increasing ε 1,2 in an oscillatory manner at even multiples of π/ω along τ 01 axis. When the delays are modulated with different modulation types, the behavior of the domain structure with increasing ε 1,2 is determined by the dominant feedback term.
Finally, one might ask the question, to which extent the extended time-delayed feedback control (ETDFC) by Socolar et al. [38,39] is also affected by a modulation on top of the discrete exponential delay distribution created by the original method. In Fig. 13 we have already presented the results for an asymmetric choice of the coupling gains K 1 = 0.05 and K 2 = 0.2. With a corresponding choice of τ 01 = 2τ 02 , which gives a special section through the shown parameter planes, the control scheme can already be seen as a rudimentary ETDFC scheme with τ 0 = τ 02 , in which the longer delays nτ 0 with n > 2 are neglected. A fully developed ETDFC scheme together with a modulation of the delays would not show any important features that are not already documented in the modulated MDFC.

VI. CONCLUSIONS
Extension of the delayed feedback control of unstable equilibria by introducing time-varying delays in one or more independent delayed feedback lines in the control circuit significantly enlarges the stability domain in the control parameter space. In addition to the domain enlargement, the variable-delay feedback method is successful in stabilizing unstable states even for those cases where the control always fails in the constant delay case.
We have shown that an analytical investigation of the control domains becomes possible in the range of highfrequency delay modulations, in which case the variabledelay term in the controller can be approximated by an equivalent distributed delay term. Our approach is motivated by both simulations and real experiments [45], which suggest that when the frequency of the delay modulation is comparable to the system frequencies, the system dynamics can be very well approximated with the distributed-delay effect. Restricting our attention to this high-frequency limit, we have investigated the control efficiency of the variable-delay feedback method by considering a simple normal form of an unstable focus and two different realizations of the controller which are experimentally relevant: a non-diagonal coupling of the control force realized via a phase-dependent coupling matrix, and the incorporation of an additional small constant delay term in all the arguments of the feedback force that represents the latency of the control-line realization. In addition, we have explored a simple realization of a multiple delay feedback controller consisting of two independent delay lines of Pyragas type with time-varying delays. In each case, the variable-delay feedback control with a finite modulation of the delay is shown superior with respect to the constant delay case. This renders the proposed control method promising for further practical implementations in real experiments.
The two-dimensional linear system used in the analysis is generic for all systems with unstable fixed points of a focus type, preserving the essential features of the higherdimensional dynamics around the equilibrium. Thus, the results obtained in this paper aim to give a succinct explanation of the delayed feedback control mechanism with time-dependent delay in real systems (regular or chaotic) with similar bifurcation properties. As in the original time-delayed feedback control, the torsion of the orbit is a necessary condition for the control method to be able to stabilize equilibria.
Since in practice the equivalence of the distributed delay case holds already for fairly low modulation frequencies, we propose that variable-delay feedback is a convenient experimental method for realizing distributed delay feedback with different delay distribution.