Dissipative systems with nonlocal delayed feedback control

Josua Grawitter; Reinier van Buel; Christian Schaaf; Holger Stark

doi:10.1088/1367-2630/aae998

1. 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 [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.

As control strategies several time-delayed feedback schemes with one or multiple delay times have been proposed [11–13]. Among these, the so-called Pyragas scheme is particularly well suited for chaotic systems and to stabilize unstable periodic orbits [14]. It has since been applied to various dynamical systems, such as lasers [15–17] and neural networks [18]. 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, the Pyragas scheme is called a noninvasive control scheme, as its stabilizing force ideally vanishes in a stabilized state [19].

Earlier theoretical studies [8, 10, 20] and experiments [9, 21] applied delayed feedback to the spatiotemporal dynamics of specific systems. One study showed that a system with two delay times can be mapped to the complex Ginzburg–Landau equation [13]. While several investigations of purely temporal systems with an intrinsic latency have focused on stabilizing unstable foci [22–24], 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 figure 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 investigate the response for conditions 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 to destabilize 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 limit the strength of the feedback control force with a sigmoid function. As a result, stable limit cycles appear in the otherwise 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 [8]. This suggests that our model captures the essential features of a spatially extended dissipative systems when subjected to nonlocal delayed feedback.

Figure 1. Refer to the following caption and surrounding text. — **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 ${\tau }_{{\rm{s}}}$ to location A.
Download figure:
Standard image High-resolution image

**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 ${\tau }_{{\rm{s}}}$ to location A.
Download figure:
Standard image High-resolution image

We present the linear dissipative model in section 2, derive the characteristic function of its fixed point, and describe our numerical methods. In section 3 we investigate and discuss the linear stability of the fixed point. In section 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 section 5.

2. Linear model with nonlocal delayed feedback

2.1. 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

$\begin{eqnarray}&&A(t)={\int }_{-\infty }^{\infty }\chi (t-{t}^{{\prime} })F({t}^{{\prime} })\,{\rm{d}}{t}^{{\prime} },\end{eqnarray} \tag{ 1 }$

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 ${\tau }_{{\rm{s}}}$ , before A is affected by F. We describe this causal link in χ using the Heaviside function ${\rm{\Theta }}(t-{\tau }_{{\rm{s}}})$ . Furthermore, in dissipative systems the response to some perturbation often decays exponentially with rate α and the response function becomes

$\begin{eqnarray}&&\chi (t)={\rm{\Theta }}(t-{\tau }_{{\rm{s}}})\exp [-\alpha (t-{\tau }_{{\rm{s}}})].\end{eqnarray} \tag{ 2 }$

Substituting χ(t) into equation (1) and differentiating both sides with respect to t, we find

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}A(t)}{{\rm{d}}t}=-\alpha A(t)+F(t-{\tau }_{{\rm{s}}}).\end{eqnarray} \tag{ 3 }$

Following Pyragas [14], we now implement for F(t) delayed feedback control with delay time ${\tau }_{{\rm{c}}}$ , control strength $\kappa$ , and a constant external force F₀,

$\begin{eqnarray}&&F(t)={F}_{0}-\kappa [A(t)-A(t-{\tau }_{{\rm{c}}})],\end{eqnarray} \tag{ 4 }$

and substitute the expression into equation (3):

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}A(t)}{{\rm{d}}t}=-\alpha A(t)+{F}_{0}-\kappa [A(t-{\tau }_{{\rm{s}}})-A(t-{\tau }_{{\rm{s}}}-{\tau }_{{\rm{c}}})].\end{eqnarray} \tag{ 5 }$

This equation has one fixed point A^*, where ${\rm{d}}A/{\rm{d}}t=0$ , which is determined by the first two terms because the delayed control term vanishes for constant solutions

$\begin{eqnarray}&&{A}^{* }=\displaystyle \frac{{F}_{0}}{\alpha }.\end{eqnarray} \tag{ 6 }$

We nondimensionalize time t and $\kappa$ with α, force F₀ with $\alpha {A}^{* }$ , and A with A^* to obtain

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}A(t)}{{\rm{d}}t}=-A(t)+1-\kappa [A(t-{\tau }_{{\rm{s}}})-A(t-{\tau }_{{\rm{s}}}-{\tau }_{{\rm{c}}})].\end{eqnarray} \tag{ 7 }$

In the following we study this form of the delay differential equation (DDE).

Note that for $\kappa \lt 0$ the delayed feedback term acts to destabilize A, because when $A(t)\gt A(t-{\tau }_{{\rm{c}}})$ the feedback term gives a positive contribution to ${\rm{d}}A/{\rm{d}}t$ . Already for ${\tau }_{{\rm{s}}}=0$ , this destabilizes the system for sufficiently large $| \kappa |$ . In the following we take $\kappa \geqslant 0$ to explore the role of ${\tau }_{{\rm{s}}}$ in the model.

2.2. Characteristic function

Within control theory, equation (7) is categorized as a closed loop system with delayed feedback. Because such systems can become unstable [25], we investigate the linear stability of the fixed point at A = 1 by introducing the perturbation ansatz

$\begin{eqnarray}&&A(t)=1+\varepsilon \exp (\lambda t)\end{eqnarray} \tag{ 8 }$

with $| \varepsilon | \ll 1$ and complex eigenvalues $\lambda \in {\mathbb{C}}$ into equation (7), which gives a transcendental equation for λ,

$\begin{eqnarray}&&\lambda =-1-\kappa {{\rm{e}}}^{-\lambda {\tau }_{{\rm{s}}}}(1-{{\rm{e}}}^{-\lambda {\tau }_{{\rm{c}}}}).\end{eqnarray} \tag{ 9 }$

In the following we solve this equation numerically and study the properties of its solutions.

The eigenvalues λ are roots of the characteristic function $g(\lambda )$ ,

$\begin{eqnarray}&&g(\lambda )=\lambda +1+\kappa {{\rm{e}}}^{-\lambda {\tau }_{{\rm{s}}}}(1-{{\rm{e}}}^{-\lambda {\tau }_{{\rm{c}}}}).\end{eqnarray} \tag{ 10 }$

Note that they cannot be expressed using the Lambert W function, as is usually possible for time-DDEs [19], 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 $\bar{\lambda }$ is a root of $g(\lambda )$ since $\kappa$ , ${\tau }_{{\rm{s}}}$ , and ${\tau }_{{\rm{c}}}$ are real parameters. Second, in appendix B we show that any root λ of g with $\mathrm{Re}(\lambda )\gt 0$ has a nonzero imaginary part bounded by $2\kappa$ , $| \mathrm{Im}(\lambda )| \lt 2\kappa$ . 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.

2.3. Numerical methods

2.3.1. Root finding

In general, the roots of the characteristic function $g(\lambda )$ 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 that the dominant eigenvalue for any parameter combination is contained in the compact complex region ${ \mathcal B }\subset {\mathbb{C}}$

$\begin{eqnarray}&&\begin{array}{l}{ \mathcal B }=\left\{\lambda \in {\mathbb{C}}| -1\leqslant \mathrm{Re}(\lambda )\leqslant \max \left[0,\displaystyle \frac{1}{{\tau }_{{\rm{s}}}}W(2\kappa {\tau }_{{\rm{s}}}{{\rm{e}}}^{{\tau }_{{\rm{s}}}})-1\right]\right.,\\ \quad 0\leqslant \mathrm{Im}(\lambda )\leqslant \kappa \exp [-{\tau }_{{\rm{s}}}\mathrm{Re}(\lambda )]\{1+\exp [-{\tau }_{{\rm{c}}}\mathrm{Re}(\lambda )]\}\}.\end{array}\end{eqnarray} \tag{ 11 }$

We search a rectangular superset of ${ \mathcal B }$ (see appendix B) using an interval Newton method [26], which provably finds all complex roots of $g(\lambda )$ in ${ \mathcal B }$ and is implemented in [27].

2.3.2. Numerical continuation

To track individual eigenvalues through parameter space, we implement a simple numerical continuation scheme for the eigenvalue equation $g(\lambda )=0$ . We explain it for the general problem of finding the roots ${\boldsymbol{u}}(p)$ of a function ${\boldsymbol{f}}({\boldsymbol{u}},p)$ . Given a set of starting values $({{\boldsymbol{u}}}_{0},{p}_{0})$ with ${\boldsymbol{f}}({{\boldsymbol{u}}}_{0},{p}_{0})=0$ , the continuation scheme first calculates an approximate step ${\tilde{{\boldsymbol{u}}}}_{1}$ along the tangent vector of the implicit function ${\boldsymbol{u}}(p)$ at ${{\boldsymbol{u}}}_{0}$ .

$\begin{eqnarray}&&{\tilde{{\boldsymbol{u}}}}_{1}={{\boldsymbol{u}}}_{0}-{\rm{\Delta }}p{\left.{({{\rm{\nabla }}}_{{\boldsymbol{u}}}{\boldsymbol{f}})}^{-1}\displaystyle \frac{\partial {\boldsymbol{f}}}{\partial p}\right|}_{({{\boldsymbol{u}}}_{0},{p}_{0})}.\end{eqnarray} \tag{ 12 }$

Here, ${\rm{\Delta }}p$ is a small number (in our calculations 10⁻²) and ${{\rm{\nabla }}}_{{\boldsymbol{u}}}{\boldsymbol{f}}$ is the Jacobian matrix of ${\boldsymbol{f}}$ w.r.t. ${\boldsymbol{u}}$ . The approximate value is then refined iteratively by solving ${\boldsymbol{f}}({{\boldsymbol{u}}}_{1},{p}_{1})=0$ for ${{\boldsymbol{u}}}_{1}$ using Newton’s method with ${p}_{1}={p}_{0}+{\rm{\Delta }}p$ and starting from ${\tilde{{\boldsymbol{u}}}}_{1}$ . The procedure produces a new set of values $({{\boldsymbol{u}}}_{1},{p}_{1})$ and is repeated as necessary to find approximate curves in parameter space. The derivatives are calculated using the library of [28] and Newton’s method is implemented in [29].

2.3.3. 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 [30] and the 4th-order Rosenbrock method RODAS [31] when using the ramp function to bound the control force (see section 4). Both methods are implemented in the software package DifferentialEquations.jl [32].

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 equation (18) in section 4) the method initially requires a history function describing A(t) for all times $-({\tau }_{{\rm{s}}}+{\tau }_{{\rm{c}}})\lt t\leqslant 0$ . The history function is used to integrate over the interval $0\lt t\leqslant {\tau }_{{\rm{s}}}$ because $A(t-{\tau }_{{\rm{s}}})$ [and $A(t-{\tau }_{{\rm{s}}}-{\tau }_{{\rm{c}}})$ ] are then known until $t={\tau }_{{\rm{s}}}$ . For remaining intervals $n{\tau }_{{\rm{s}}}\lt t\leqslant (n+1){\tau }_{{\rm{s}}}$ (with $n\in {\mathbb{N}}$ ) integration continues using A(t) from the previous intervals.

2.3.4. Periodic solutions

To find time-periodic solutions with period T for the DDE

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}A(t)}{{\rm{d}}t}=f[A(t),A(t-{\tau }_{1}),A(t-{\tau }_{2})]\end{eqnarray} \tag{ 13 }$

with delays ${\tau }_{1}$ and ${\tau }_{2}$ for a given function f and initial conditions, we use the following method. First, we construct a truncated Fourier series $\tilde{A}(t)$ with N real coefficients corresponding to the N lowest-order sine and cosine harmonics with periodicity T. Initially, the coefficients are chosen such that $\tilde{A}(t)=0.25\sin (2\pi t/T)$ . From the given DDE we introduce a residual function $\varepsilon (t)$ , which describes the error of approximation $\tilde{A}$ ,

$\begin{eqnarray}&&\varepsilon (t)=\displaystyle \frac{{\rm{d}}\tilde{A}(t)}{{\rm{d}}t}-f[\tilde{A}(t),\tilde{A}(t-{\tau }_{1}),\tilde{A}(t-{\tau }_{2})].\end{eqnarray} \tag{ 14 }$

We also truncate $\varepsilon (t)$ in Fourier space with N real coefficients (corresponding to the same modes as before). We apply a root-finding method to find the set of coefficients for $\tilde{A}(t)$ , which minimizes the absolute values of the coefficients of $\varepsilon (t)$ . When all coefficients of $\varepsilon (t)$ become vanishingly small, we have arrived at a good approximation for A(t). To discretize both $\tilde{A}$ and ε in Fourier space, we use the library of [33] with N = 20. For root finding we use up to 10 iterations of a trust region method implemented in [29].

Note that for many periods T the only periodic solutions are fixed points. To distinguish fixed points from limit cycles, we use the ratio of the amplitudes of $\tilde{A}(t)$ and ε(t): because fixed points have exactly zero amplitude but the residual nevertheless fluctuates, the amplitude ratio is small for fixed points and large for limit cycles. By maximizing the ratio we obtain the periods of limit cycles in our system. To perform the numerical optimization of T, we use a golden-section search implemented in [34].

Finally, we test the stability of limit cycles numerically using the time-integration methods described in section 2.3.3 by introducing them as history functions and initial conditions. Due to numerical error, unstable limit cycles eventually diverge from their orbits while stable limit cycles remain there. We observe trajectories for 100 units of time to determine stability.

All numerical calculations are performed using the Julia programming language [35] and all plots are created using the matplotlib package [36].

3. Linear stability analysis

3.1. Stability for vanishing delays

In the limiting case of vanishing system delay, ${\tau }_{{\rm{s}}}\to 0$ , the fixed point is always stable. To prove this, we set ${\tau }_{{\rm{s}}}=0$ and take the real part of g:

$\begin{eqnarray}&&\mathrm{Re}[g(\lambda )]=\mathrm{Re}(\lambda )+1+\kappa -\kappa {{\rm{e}}}^{-{\tau }_{{\rm{c}}}\mathrm{Re}(\lambda )}\cos [{\tau }_{{\rm{c}}}\mathrm{Im}(\lambda )].\end{eqnarray} \tag{ 15 }$

We prove by contradiction¹ : suppose that $\mathrm{Re}(\lambda )\geqslant 0$ together with $\kappa ,{\tau }_{{\rm{c}}}\gt 0$ . Then the final term in equation (15) is the only (potentially) negative contribution necessary to obtain $\mathrm{Re}(g)=0$ . However, because the absolute value of this term is always smaller than or equal to $\kappa$ , roots with $\mathrm{Re}(\lambda )\geqslant 0$ cannot exist. Therefore, all eigenvalues have $\mathrm{Re}(\lambda )\lt 0$ and the fixed point must be stable for ${\tau }_{{\rm{s}}}\to 0$ .

In the limit of vanishing control delay, ${\tau }_{{\rm{c}}}\to 0$ , the fixed point is stable because the characteristic function

$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{\tau }_{{\rm{c}}}\to 0}g(\lambda )=\lambda +1\end{eqnarray} \tag{ 16 }$

only has one (real) root at $\lambda =-1$ . This is also clear from equation (7) since the delay term vanishes completely.

3.2. Stability-to-instability transition

A fixed point is unstable if $\mathrm{Re}(\lambda )\gt 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 ${\tau }_{{\rm{s}}}\gt 0$ and ${\tau }_{{\rm{c}}}\gt 0$ , there is a critical control strength ${\kappa }_{\mathrm{crit}}({\tau }_{{\rm{s}}},{\tau }_{{\rm{c}}})$ determined by $\mathrm{Re}(\lambda )=0$ for the dominant eigenvalue and $\mathrm{Re}(\lambda )\leqslant 0$ for all others beyond which the fixed point is unstable for all $\kappa \gt {\kappa }_{\mathrm{crit}}$ . The ${\kappa }_{\mathrm{crit}}$ form a manifold in parameter space; cross sections for various ${\tau }_{{\rm{s}}}$ are displayed in figure 2. For ${\tau }_{{\rm{c}}}\to 0$ the critical value ${\kappa }_{\mathrm{crit}}$ diverges like ${\tau }_{{\rm{c}}}^{-1}$ . In this limit $g(\lambda )$ is approximated to linear order in ${\tau }_{{\rm{c}}}$ as

$\begin{eqnarray}&&g(\lambda )=\lambda +1+\kappa {\tau }_{{\rm{c}}}\lambda \,{{\rm{e}}}^{-\lambda {\tau }_{{\rm{s}}}},\end{eqnarray} \tag{ 17 }$

from which follows that any root λ stays the same as long as the product $\kappa {\tau }_{{\rm{c}}}$ remains constant. For increasing ${\tau }_{{\rm{c}}}$ specific values of ${\kappa }_{\mathrm{crit}}$ exist, where a different eigenvalue becomes dominant. We observe that the resulting cusps become less prominent as ${\tau }_{{\rm{c}}}$ increases and eventually ${\kappa }_{\mathrm{crit}}$ as a function of ${\tau }_{{\rm{c}}}$ approaches a constant value for each ${\tau }_{{\rm{s}}}$ . Notably, these cusps imply that for $\kappa$ close to the cusp value alternating regions of stability and instability occur as ${\tau }_{{\rm{c}}}$ increases. We also observe that for small ${\tau }_{{\rm{s}}}$ the cusps are closer to each other w.r.t. ${\tau }_{{\rm{c}}}$ and ${\kappa }_{\mathrm{crit}}$ diverges with decreasing ${\tau }_{{\rm{s}}}$ . From these observations and the findings of section 3.1, we conclude that only the combination of both nonzero delays causes the instability.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Stability–instability transition curves in the $\kappa$ – ${\tau }_{{\rm{c}}}$ plane for various ${\tau }_{{\rm{s}}}$ with the stable region to the left of each curve.
Download figure:
Standard image High-resolution image

The dominant eigenvalue displays some general features (shown in figure 3): its real part is smallest (i.e. close to −1) in regions with either $\kappa$ or ${\tau }_{{\rm{c}}}$ close to zero, where the control term becomes negligible. Conversely, this implies that delayed feedback slows down the decay of individual modes. The imaginary part of the dominant eigenvalue is nonzero and displays clear discontinuities w.r.t. ${\tau }_{{\rm{c}}}$ and $\kappa$ in the unstable region and close to the transition curve ${\kappa }_{\mathrm{crit}}$ (see figure 3(b)). Since also $\mathrm{Re}(\lambda )=0$ at the transition, there is a Hopf bifurcation. Furthermore, as ${\tau }_{{\rm{c}}}$ increases the imaginary part repeatedly decreases and then jumps to a larger value when a different eigenvalue becomes dominant. Figure 3(a) shows that the jump lines correspond to valleys in the real part of the dominant eigenvalue.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Dominant eigenvalue λ color-coded in ${\tau }_{{\rm{c}}}$ – $\kappa$ plane together with stability–instability transition curve ${\kappa }_{\mathrm{crit}}$ (dashed red/white lines) for system delay ${\tau }_{{\rm{s}}}=0.2$ : (a) real part and (b) imaginary part.
Download figure:
Standard image High-resolution image

More generally, ${\kappa }_{\mathrm{crit}}({\tau }_{{\rm{c}}})$ is the envelope of a family of Hopf curves, i.e. lines in parameter space corresponding to $\mathrm{Re}({\lambda }_{\pm })=0$ for each conjugate pair of complex eigenvalues λ_±. We determine the Hopf curves numerically as implicit functions ${\boldsymbol{u}}({\tau }_{{\rm{c}}})$ defined by $g({\boldsymbol{u}},{\tau }_{{\rm{c}}})=0$ with ${\boldsymbol{u}}=(\mathrm{Im}(\lambda ),\kappa )$ and with constant ${\tau }_{{\rm{s}}}$ (note that g is complex and therefore constrains both coordinates of ${\boldsymbol{u}}$ ). Six Hopf curves for successive eigenvalues are displayed in figure 4 with ${\tau }_{{\rm{s}}}=0.2$ (corresponding to figure 3). Extending the Hopf curves to large $\kappa$ , we find that they approach constant values for ${\tau }_{{\rm{c}}}$ . The cusps visible in figures 2 and 3 correspond to intersecting Hopf curves, i.e. two pairs of eigenvalues, where the real part switches from negative to positive at the same set of parameters. The pairs of eigenvalues have positive real parts on the right-hand side of their respective curves and, therefore, crossings imply that several eigenvalues can have positive real parts for the same set of parameters. Furthermore, when we project the numerically calculated curves $\kappa ({\tau }_{{\rm{c}}})$ onto the $\mathrm{Im}(\lambda )$ – $\kappa$ plane, they all fall onto the same line (see inset of figure 3). This suggests there is a global minimum of critical control strength ${\kappa }_{\mathrm{crit}}^{\min }$ for the envelope ${\kappa }_{\mathrm{crit}}({\tau }_{{\rm{c}}})$ below which no instability can occur. Lastly, we observe that successive Hopf curves $\kappa ({\tau }_{{\rm{c}}})$ become flatter, which means their envelope ${\kappa }_{\mathrm{crit}}({\tau }_{{\rm{c}}})$ approaches a straight line for large ${\tau }_{{\rm{c}}}$ .

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Hopf curves of six successive eigenvalues λ (solid colored lines) in the ${\tau }_{{\rm{c}}}$ – $\kappa$ plane for system delay ${\tau }_{{\rm{s}}}=0.2$ . Inset: in the $\mathrm{Im}(\lambda )$ – $\kappa$ plane all six Hopf curves project onto the same curve (solid black line). As a result, they all share the same minimum of critical control strength ${\kappa }_{\mathrm{crit}}^{\min }$ (dashed gray line).
Download figure:
Standard image High-resolution image

We close with two comments. First, we note the similarities between our dominant eigenvalues and eigenvalues described in [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 equation (10). Furthermore, the diffusive response function (compare equation (20) in appendix A with $\alpha \to 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 equation (2).

Second, when Pyragas feedback is applied to a system with spatial components like the Swift–Hohenberg equation [37, 38] or an reaction-diffusion system [39], the relation $\kappa {\tau }_{{\rm{c}}}=1$ marks the instability boundary, where spontaneous motion of spatial patterns occurs. In our simple model we do not find an explicit expression for the transition line ${\kappa }_{\mathrm{crit}}({\tau }_{{\rm{c}}})$ but a similar scaling in case of small control delays ${\tau }_{{\rm{c}}}\to 0$ (see above).

4. Bounded control force

Typically, the response of a physical system to an external stimulus is bounded by nonlinearities because ultimately no control mechanism can exert an infinite force. 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 [40], or for hydrodynamic vortex diffusion in a circular geometry to stabilize unstable modes [8].

4.1. Model and setup

We start by modifying the DDE of equation (7) such that the time-delayed control term is limited by a monotonically increasing odd function $\sigma (x)$ with $\sigma (\pm \infty )=\pm 1$ and $\sigma ^{\prime} (0)=1$ :

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}A(t)}{{\rm{d}}t}=-A(t)-\sigma \{\kappa [A(t-{\tau }_{{\rm{s}}})-A(t-{\tau }_{{\rm{s}}}-{\tau }_{{\rm{c}}})]\}.\end{eqnarray} \tag{ 18 }$

Compared to equation (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 equation (7) without the term +1.

To realize upper and lower bounds for the control forces, we consider three sigmoid functions: the hyperbolic tangent $\sigma (x)=\tanh (x)$ , which approaches ±1 exponentially, an algebraic sigmoid $\sigma (x)=x{(1+| x| )}^{-1}$ , which approaches ±1 like a power law, and the nonsmooth ramp function $\sigma (x)=\mathrm{ramp}(x)=\min [1,\max (-1,x)]$ . All three are displayed in figure 5.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Three realizations of a sigmoid function for implementing delayed control forces with upper and lower bounds plotted versus x. Three asymptotic scalings for large $| x|$ are realized: flat (ramp), exponential (tanh), and power law.
Download figure:
Standard image High-resolution image

$| x| $ — **Figure 5.** Three realizations of a sigmoid function for implementing delayed control forces with upper and lower bounds plotted versus x. Three asymptotic scalings for large $| x|$ are realized: flat (ramp), exponential (tanh), and power law.
Download figure:
Standard image High-resolution image

We solve equation (18) numerically² as described in section 2.3 with the history function $A(t\lt 0)=0$ and a small initial perturbation³ , A(t = 0) = 10⁻³. 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 section 2 by calculating the time evolution of A(t) until time T = 10³ and by examining the frequencies and amplitudes of the occurring limit cycles for times 0.8T < t < T. We chose these initial conditions because we are specifically interested in the long-time behavior following a simple perturbation from the fixed point without any prescribed mode or frequency.

4.2. 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 figure 6). 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 [14], because the feedback term does not vanish when the limit cycle is reached. Generally, the limit cycles for all sigmoid functions should converge for $\kappa \to \infty$ 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 figure 6 with figure 3(b). In particular, the jumps from one mode to the other agree well for both frequencies. However, the amplitude ${\rm{\Delta }}A$ of the limit cycle ${A}_{\mathrm{lc}}(t)$ defined as ${\rm{\Delta }}A=[{\max }_{t}\{{A}_{\mathrm{lc}}(t)\}-{\min }_{t}\{{A}_{\mathrm{lc}}(t)\}]/2$ is generally not related to the real part of the dominant eigenvalue (compare amplitudes in the right column of figure 6 with figure 3(a)). It rather behaves like the frequency of the limit cycle with one difference: it increases as ${\tau }_{{\rm{c}}}$ grows and drops to a smaller value when the long-time dynamics switches to a different limit cycle. So, at the discontinuity lines the sudden increase in frequency is accompanied by a sudden decrease in amplitude. This is explicitly demonstrated in figure 7.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Frequencies and amplitudes of the limit cycles resulting from the DDE of equation (18) color-coded in the $\kappa$ – ${\tau }_{{\rm{c}}}$ parameter space at ${\tau }_{{\rm{s}}}=0.2$ for different realizations σ of the bounded control force: (a) ramp sigmoid, (b) hyperbolic tangent, and (c) algebraic sigmoid. Zero frequency (dark blue) indicates the region of stable fixed point (no oscillations). The dashed white line indicates the transition from stable to unstable fixed point.
Download figure:
Standard image High-resolution image

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Limit cycles of a system with feedback control bounded by the hyperbolic tangent as sigmoid function at two different control delays ${\tau }_{{\rm{c}}}$ for $\kappa =5$ and ${\tau }_{{\rm{s}}}=0.2$ . The system with larger ${\tau }_{{\rm{c}}}$ has larger periodicity but smaller amplitude.
Download figure:
Standard image High-resolution image

Using the method described in section 2.3.4 we analyzed the two limit cycles displayed in figure 7 in more detail. We find that they co-exist for both sets of parameters (with slightly shifted amplitudes and frequencies). For ${\tau }_{{\rm{c}}}=1.0$ only the limit cycle with the larger frequency is stable, while for ${\tau }_{{\rm{c}}}=1.2$ both limit cycles are stable and the long-time dynamics depends on initial conditions. Apparently, in the bistable case the high-frequency limit cycle attracts trajectories departing from the unstable fixed point, which where attracted to the low-frequency limit cycle for ${\tau }_{{\rm{c}}}=1.0$ . Notably, the transition to bistability occurs close to one of the Hopf curves displayed in figure 4, which suggests that the second limit cycle becomes stable as an additional pair of complex eigenvalues acquires positive real parts. We, therefore, expect that for larger ${\tau }_{{\rm{c}}}$ even more than two stable limit cycles co-exist, depending on the number of complex eigenvalues with positive real part⁴ .

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 figure 6) 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 $\kappa$ . 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 figure 3(b). However, there are slight deviations for the smooth sigmoid functions.

4.3. Transient pulse trains

In particular, for large control times ${\tau }_{{\rm{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 figure 8(a). These pulse trains grow into regular limit cycles, as shown in figure 8(b), while for stable fixed points their amplitude decays to zero. Generally, we observe that pulse trains repeat with a periodicity given by ${\tau }_{{\rm{c}}}$ and their oscillation frequency is close to the imaginary part of the dominant eigenvalue $\mathrm{Im}(\lambda )/2\pi$ .

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Dynamic response of the system to feedback control bounded by the hyperbolic tangent function for ${\tau }_{{\rm{s}}}=0.4$ , ${\tau }_{{\rm{c}}}=25$ , and $\kappa =2.1$ : (a) transient pulse trains and (b) limit cycle oscillations (note the different scales).
Download figure:
Standard image High-resolution image

**Figure 8.** Dynamic response of the system to feedback control bounded by the hyperbolic tangent function for ${\tau }_{{\rm{s}}}=0.4$ , ${\tau }_{{\rm{c}}}=25$ , and $\kappa =2.1$ : (a) transient pulse trains and (b) limit cycle oscillations (note the different scales).
Download figure:
Standard image High-resolution image

5. 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 [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 stability–instability 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 [14], 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, which attracts trajectories starting from the fixed point, match well the imaginary part of the fixed point’s dominant eigenvalue. Discontinuity lines are visible, which occur when the long-time behavior switches to a different limit cycle causing frequency to jump up and amplitude to drop 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 nontrivial behavior. It also provides an example for double delay systems, which have recently attracted attention [42–44]. In future studies we will apply delayed feedback to specific nonlinear soft matter systems such as photoresponsive fluid interfaces [45] and viscoelastic flow in Taylor–Couette geometry [46]. The work presented in this article will help us to categorize the observed spatiotemporal dynamics.

Acknowledgments

We acknowledge financial support from DFG (German Research Foundation) via International Research Training Group 1524 and Collaborative Research Center 910. We further acknowledge support from the DFG and the Open Access Publication Funds of TU Berlin.

Appendix A.: A physical system with intrinsic delay

We consider a diffusion-reaction equation of a substance with density $\rho ({\boldsymbol{r}},t)$ , which decays at a rate α in two-dimensions:

$\begin{eqnarray}&&{\partial }_{t}\rho ({\boldsymbol{r}},t)=D{{\rm{\nabla }}}^{2}\rho ({\boldsymbol{r}},t)-\alpha \rho (x,t),\end{eqnarray} \tag{ 19 }$

where D is the diffusion constant. As initial condition at t = 0 we choose a delta peak located at ${\boldsymbol{r}}=0$ , $\rho ({\boldsymbol{r}},0)=N\delta ({\boldsymbol{r}})$ . The solution of this problem is given by

$\begin{eqnarray}&&\rho ({\boldsymbol{r}},t)=\displaystyle \frac{N}{4\pi {Dt}}\exp \left(-\alpha t-\displaystyle \frac{{r}^{2}}{4{Dt}}\right),\end{eqnarray} \tag{ 20 }$

where $r=| {\boldsymbol{r}}|$ is the distance from the perturbation. The density at any point ${{\boldsymbol{r}}}_{0}\ne 0$ increases up to the time

$\begin{eqnarray}&&{t}_{0}=\displaystyle \frac{1}{2\alpha }\left(\sqrt{1+\displaystyle \frac{\alpha {r}_{0}^{2}}{D}}-1\right)\end{eqnarray} \tag{ 21 }$

and then decreases to zero. Note that for pure diffusion ( $\alpha \to 0$ ) the maximum is reached at ${t}_{0}={r}_{0}^{2}{(4D)}^{-1}$ . Thus, a disturbance initiated at r = 0 needs the intrinsic delay time t₀ to reach ${{\boldsymbol{r}}}_{0}$ . To approximate this behavior in a response function of the form given in equation (2), we assume a step function where the substance ρ jumps to its maximum value $\rho ({{\boldsymbol{r}}}_{0},{t}_{0})$ and then decays exponentially

$\begin{eqnarray}&&\chi ({{\boldsymbol{r}}}_{0},t)={\rm{\Theta }}(t-{t}_{0})\rho ({{\boldsymbol{r}}}_{0},{t}_{0})\exp [-\alpha (t-{t}_{0})].\end{eqnarray} \tag{ 22 }$

Thus, ${\tau }_{{\rm{s}}}={t}_{0}$ is the intrinsic delay time and the physical decay rate α of our substance enters directly the response function. Because the step function prevents any response for $t\lt {t}_{0}$ , t₀ is an effective propagation time.

Appendix 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 $\mathrm{Re}[g(\lambda )]=0$ :

$\begin{eqnarray}&&\displaystyle \frac{\mathrm{Re}(\lambda )+1}{\kappa }=\exp [-({\tau }_{{\rm{s}}}+{\tau }_{{\rm{c}}})\mathrm{Re}(\lambda )]\cos [({\tau }_{{\rm{s}}}+{\tau }_{{\rm{c}}})\mathrm{Im}(\lambda )]-\exp [-{\tau }_{{\rm{s}}}\mathrm{Re}(\lambda )]\cos [{\tau }_{{\rm{s}}}\mathrm{Im}(\lambda )].\end{eqnarray} \tag{ 23 }$

The cosines may at most change the signs of the terms such that both contribute positively. Thus

$\begin{eqnarray}&&\displaystyle \frac{\mathrm{Re}(\lambda )+1}{\kappa }\leqslant \exp [-{\tau }_{{\rm{s}}}\mathrm{Re}(\lambda )]+\exp [-({\tau }_{{\rm{s}}}+{\tau }_{{\rm{c}}})\mathrm{Re}(\lambda )].\end{eqnarray} \tag{ 24 }$

Assuming $\mathrm{Re}(\lambda )\geqslant 0$

$\begin{eqnarray}&&\displaystyle \frac{\mathrm{Re}(\lambda )+1}{\kappa }\leqslant 2\exp [-{\tau }_{{\rm{s}}}\mathrm{Re}(\lambda )],\end{eqnarray} \tag{ 25 }$

and we solve for $\mathrm{Re}(\lambda )$ using Lambert’s W function

$\begin{eqnarray}&&\mathrm{Re}(\lambda )\leqslant \displaystyle \frac{1}{{\tau }_{{\rm{s}}}}W(2\kappa {\tau }_{{\rm{s}}}{{\rm{e}}}^{{\tau }_{{\rm{s}}}})-1.\end{eqnarray} \tag{ 26 }$

With our previous assumption we have

$\begin{eqnarray}&&\mathrm{Re}(\lambda )\leqslant \max \left[0,\displaystyle \frac{1}{{\tau }_{{\rm{s}}}}W(2\kappa {\tau }_{{\rm{s}}}{{\rm{e}}}^{{\tau }_{{\rm{s}}}})-1\right].\end{eqnarray} \tag{ 27 }$

Notably the upper bound is independent of ${\tau }_{{\rm{c}}}$ .

Second, we find a lower bound for $\mathrm{Re}(\lambda )$ such that our search region contains at least one eigenvalue. To simplify our search, we concentrate on $\lambda \in {\mathbb{R}}$ , which is determined by

$\begin{eqnarray}&&\lambda =-1-\exp (-{\tau }_{{\rm{s}}}\lambda )[1-\exp (-{\tau }_{{\rm{c}}}\lambda )].\end{eqnarray} \tag{ 28 }$

On the interval $(-\infty ,0]$ the lhs of this equation goes continuously from $-\infty$ to 0 and the rhs goes continuously from $+\infty$ to −1. Because both sides are continuous and strictly monotonic for ${\tau }_{{\rm{c}}},{\tau }_{{\rm{s}}}\gt 0$ , they share exactly one value in the overlap, i.e. $-1\leqslant \lambda \leqslant 0$ . Therefore, a search region with $\mathrm{Re}(\lambda )\geqslant -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 $\mathrm{Im}[g(\lambda )]=0$

$\begin{eqnarray}&&\mathrm{Im}(\lambda )=\kappa \exp [-{\tau }_{{\rm{s}}}\mathrm{Re}(\lambda )]\sin [{\tau }_{{\rm{s}}}\mathrm{Im}(\lambda )]-\kappa \exp [-({\tau }_{{\rm{s}}}+{\tau }_{{\rm{c}}})\mathrm{Re}(\lambda )]\sin [({\tau }_{{\rm{s}}}+{\tau }_{{\rm{c}}})\mathrm{Im}(\lambda )].\end{eqnarray} \tag{ 29 }$

Here, we simply drop the sines, assuming all terms contribute positively

$\begin{eqnarray}&&\mathrm{Im}(\lambda )\leqslant \kappa \exp [-{\tau }_{{\rm{s}}}\mathrm{Re}(\lambda )]\{1+\exp [-{\tau }_{{\rm{c}}}\mathrm{Re}(\lambda )]\}.\end{eqnarray} \tag{ 30 }$

The rhs depends on $\mathrm{Re}(\lambda )$ and on the interval determined for $\mathrm{Re}(\lambda )$ . It is maximal for $\mathrm{Re}(\lambda )=-1$ implying $\mathrm{Im}(\lambda )\leqslant \kappa \,{{\rm{e}}}^{{\tau }_{{\rm{s}}}}[1+{{\rm{e}}}^{{\tau }_{{\rm{c}}}}]$ . Note that the symmetry $g(\bar{\lambda })=\bar{g}(\lambda )$ implies that we only need to search the positive half of the complex plane, i.e. $\mathrm{Im}(\lambda )\geqslant 0$ , because all complex roots appear in conjugate pairs.

In summary, we have shown that the following set of complex numbers ${ \mathcal B }$ must always contain the eigenvalue with the largest real part:

$\begin{eqnarray}\begin{array}{rcl}{ \mathcal B }=\{\lambda \in {\mathbb{C}}| -1 & \leqslant & \mathrm{Re}(\lambda )\leqslant \max \left[0,\displaystyle \frac{1}{{\tau }_{{\rm{s}}}}W(2\kappa {\tau }_{{\rm{s}}}{{\rm{e}}}^{{\tau }_{{\rm{s}}}})-1\right],\\ 0 & \leqslant & \mathrm{Im}(\lambda )\leqslant \kappa \exp [-{\tau }_{{\rm{s}}}\mathrm{Re}(\lambda )]\{1+\exp [-{\tau }_{{\rm{c}}}\mathrm{Re}(\lambda )]\}\}.\end{array}\end{eqnarray} \tag{ 31 }$

For simplicity and safety we search a bounded rectangular region which is a superset of ${ \mathcal B }$ :

$\begin{eqnarray}\begin{array}{rcl}{ \mathcal B }\subseteq \{\lambda \in {\mathbb{C}}| -1.1 & \leqslant & \mathrm{Re}(\lambda )\leqslant \max \left[0,\displaystyle \frac{1.1}{{\tau }_{{\rm{s}}}}W(2\kappa {\tau }_{{\rm{s}}}{{\rm{e}}}^{{\tau }_{{\rm{s}}}})-1.1\right],\\ -0.1 & \leqslant & \mathrm{Im}(\lambda )\leqslant 1.1\,\kappa {{\rm{e}}}^{{\tau }_{{\rm{s}}}}(1+{{\rm{e}}}^{{\tau }_{{\rm{c}}}})\}.\end{array}\end{eqnarray} \tag{ 32 }$

Both regions are displayed in figure B1.

Figure B1. Refer to the following caption and surrounding text. — **Figure B1.** Example of a region, which with certainty contains the dominant eigenvalue as determined in equation (31) (red, transparent), numerical search rectangle from equation (32) (white, transparent), and numerically located eigenvalues (black crosses) for $\kappa =5$ , ${\tau }_{{\rm{c}}}=1$ , ${\tau }_{{\rm{s}}}=0.4$ . The upper bound for the real parts of eigenvalues with large imaginary parts is displayed as a white, dashed line. The background shading indicates the absolute value of the characteristic function $| g(\lambda )|$ .
Download figure:
Standard image High-resolution image

**Figure B1.** Example of a region, which with certainty contains the dominant eigenvalue as determined in equation (31) (red, transparent), numerical search rectangle from equation (32) (white, transparent), and numerically located eigenvalues (black crosses) for $\kappa =5$ , ${\tau }_{{\rm{c}}}=1$ , ${\tau }_{{\rm{s}}}=0.4$ . The upper bound for the real parts of eigenvalues with large imaginary parts is displayed as a white, dashed line. The background shading indicates the absolute value of the characteristic function $| g(\lambda )|$ .
Download figure:
Standard image High-resolution image

Note that the shape of ${ \mathcal B }$ immediately implies that any eigenvalue λ with $\mathrm{Re}(\lambda )\gt 0$ has $| \mathrm{Im}(\lambda )| \lt 2\kappa$ . Furthermore, it can be shown that $\mathrm{Re}(\lambda )\gt 0$ implies $| \mathrm{Im}(\lambda )| \ne 0$ because equation (28) has no solutions for $\lambda \gt 0$ . In this case the lhs of equation (28) is always positive while the rhs is always smaller than −1 and therefore no solution with $\mathrm{Re}(\lambda )\gt 0$ and $| \mathrm{Im}(\lambda )| =0$ exists.