Complete characterization of the asymptotic stability region for single-delay second-order linear systems

1 Introduction

Consider the problem of stabilizing a second-order plant under mixed local and delayed feedback, a situation that arises in vibration suppression, machine-tool chatter, human balance control, and networked control, as presented in Equation 1:

x ̈ t + c x ̇ t + ω 0 2 x t + k d x ̇ t + k p x t + k ̄ d x ̇ t − 1 + k ̄ p x t − 1 = 0 . (1)

The corresponding characteristic equation is

s 2 + c + k d s + ω 0 2 + k p + k ̄ d s + k ̄ p e − s = 0 .

For fixed delayed-channel parameters ( a 01 , a 11 ) = ( k ̄ p , k ̄ d ), our main result yields a complete gain chart in the ( a 00 , a 10 )-plane: it characterizes the full region of stabilizing local parameters ( a 00 , a 10 ) = ( ω 0 2 + k p , c + k d ) and its boundary, rather than testing stability at isolated points. The dual “delay-margin” viewpoint—admissible τ for fixed gains—is recovered by the scaling map in Remark 11: varying τ traces a ray in parameter space, and stability holds for those τ such that the scaled point ( a 00 , a 10 ) = ( τ 2 a ̄ 00 , τ a ̄ 10 ) lies in the region D associated with ( a 01 , a 11 ) = ( τ 2 a ̄ 01 , τ a ̄ 11 ).

Despite the classical status of this problem and the extensive literature on delay-system stability—with foundational contributions by Pontryagin (1955), Ernest Bellman and Kenneth (1963), and Hale and Verduyn Lunel (1993), along with the monograph of Stépán (1989) on retarded mechanical systems—a complete, necessary, and sufficient characterization of the stabilizing set in the full four-parameter space of Equation 1 has, to our knowledge, not been reported. The difficulty is that even a single delay renders the characteristic equation transcendental, yielding infinitely many roots whose distribution must be controlled simultaneously as four coefficients vary.

A large body of work has developed computational frameworks that address this difficulty in the dual parametrization. The τ-decomposition method—rooted in the D-partition technique by Neĭmark (1949) and developed for delay equations by Cooke and Grossman (1982) and Kenneth (1986)—together with the direct method of Walton and Marshall (1987), the cluster treatment of characteristic roots (Olgac and Sipahi, 2002; Sipahi et al., 2011), and frequency-sweeping tests (Chen et al., 1995), all determine, for fixed coefficients, the set of admissible delays or crossing frequencies; comprehensive expositions are provided by Gu et al. (2003) and Michiels and Niculescu (2014). The dual problem addressed here—characterizing the stabilizing set of coefficients for a fixed delay—has received considerably less attention.

The field has remained very active over the last decade. On the computational side, pseudospectral discretization of the solution operator or the infinitesimal generator now provides efficient approximations of the rightmost characteristic roots (Breda et al., 2015), complementing large-scale quasi-polynomial root mapping (Vyhlídal and Zítek, 2009) and the eigenvalue-based framework of Michiels and Niculescu (2014). In the time domain, necessary and sufficient stability tests expressed through the delay Lyapunov matrix have reached maturity; see the recent survey (Mondié et al., 2022). Frequency-domain criteria continue to be refined from new perspectives on the characteristic equation (Tessina, 2023), and the geometric stability-switch methodology of Cooke–van den Driessche type has been extended to two delays and delay-dependent coefficients (Qi et al., 2019). A distinct recent line of research exploits multiplicity-induced dominancy (MID)—the fact that characteristic roots of maximal multiplicity of low-order retarded equations are dominant—to assign the spectral abscissa directly (Mazanti et al., 2021). All of these approaches, such as the abovementioned classical approaches, operate primarily in the dual parametrization (fixed coefficients, varying delay, or spectrum); the problem of describing exactly the stabilizing set of coefficients has remained largely open.

For second-order single-delay equations specifically, the scalar (first-order) case is classical, but the second-order case has remained open in full generality. Cahlon and Schmidt (2004) derived sufficient stability criteria using Pontryagin’s theory of quasi-polynomials under the sign-product restrictions a 00 a 01 < 0 and a 10 a 11 ≥ 0; they subsequently extended their analysis to delay differential equations of arbitrary even and odd order (Cahlon and Schmidt, 2007a; Cahlon and Schmidt, 2007b). Brethé (1994) addressed stability and stabilizability for several second-order sub-cases, while the classical problem of stabilizing an oscillator via delayed position feedback was treated by Abdallah et al. (1993). In the PID-stabilization literature, Silva et al. (2004) characterized the set of stabilizing gains for first-order and certain second-order plants using generalized Hermite–Biehler arguments. More recently, the delay margin achievable by Proportional-Integral-Derivative (PID) control of low-order plants was characterized by Ma and Chen (2019), the exact delay margin of PID control of second-order systems was obtained by Chen et al. (2022) via a pseudoconcavity argument, and strongly stabilizing PID controllers for general delay systems were designed by Appeltans et al. (2022); these studies optimize a margin along the delay axis for a structured controller, rather than describing the stabilizing region in the coefficient space of Equation 3. In the machining dynamics literature, the classical stability lobe diagrams (Altintaş and Budak, 1995; Tamás Insperger and Gábor Stépán, 2011) provide a graphical stabilizing set for the undamped sub-case a 10 = a 11 = 0. All of these results either treat restricted parameter sub-classes or provide only sufficient conditions. A complete, necessary, and sufficient characterization of the stability region in the parametric space for the full four-parameter class of second-order single-delay Equation 3 closes this gap and unifies the cited special cases (see Table 1).

TABLE 1. Summary of available stability criteria related to Equation 3. The first block collects results for the full model under additional parameter restrictions; the second block lists results for true subcases obtained by coefficient cancellations.

1.1 Contribution

In contrast to the numerical methods mentioned above, which compute stability boundaries by tracking root crossings or discretizing the solution operator, the present study provides an analytical and constructive characterization of the entire asymptotic stability region in the parameter space for the general second-order single-delay linear system (Equation 2):where all coefficients are real.

x ̈ t + a ̄ 10 x ̇ t + a ̄ 00 x t + a ̄ 11 x ̇ t − τ + a ̄ 01 x t − τ = 0 , τ ∈ R > 0 , (2)

Using the standard time-scaling normalization to unit delay, we reduce the analysis to the system shown in Equation 3:which depends on four real parameters a 00 , a 10 , a 01 , a 11. We then characterize the set of all parameter quadruples ( a 00 , a 10 , a 01 , a 11 ) for which the trivial solution of Equation 3 is asymptotically stable.

x ̈ t + a 10 x ̇ t + a 00 x t + a 11 x ̇ t − 1 + a 01 x t − 1 = 0 , (3)

Specifically, for fixed values of ( a 01 , a 11 ), we define a region D in the ( a 00 , a 10 )-plane via the supremum envelope (Definition 2) and prove that the zero solution of Equation 3 is asymptotically stable if and only if ( a 00 , a 10 ) ∈ D (Theorem 29). In contrast to classical “delay-margin” analyses that focus on admissible intervals of τ for a fixed set of coefficients, our result directly characterizes an entire stabilizing region in coefficient space; delay margins and gain–delay charts can then be recovered by the scaling map in Remark 11. This characterization

• unifies and extends several classical stability criteria for special cases of Equation 3 (summarized in Table 1), providing necessary and sufficient conditions without auxiliary parameter restrictions;
• enables robustness margins defined as distances to the boundary of D (or, equivalently, to Γ);
• facilitates synthesis and tuning: D is an epigraph, and its boundary is the computable curve a 10 = A ̂ 10 ( a 00 );
• can yield simple closed-form conditions for structured parameter relationships, as illustrated in our examples;
• is directly computable: Algorithm 1 provides a numerical recipe to evaluate A ̂ 10 ( a 00 ) and to sample the boundary of D, enabling standard visualization and design workflows (including delay-margin computations via Remark 11).

Moreover, for fixed ( a 01 , a 11 ), we show that D is an epigraph, hence connected and straightforward to visualize.

It is worth delineating precisely which ingredients of this characterization are classical and which are new. The parametrization of the imaginary-axis-crossing locus (Lemma 2), the invariance of the number of unstable roots on each connected component of the complement of the crossing curve (Corollary 14), and the transversal crossing-direction computations (Corollaries 7, 8) are standard elements of the D-partition method (Neĭmark, 1949; Cooke and Grossman, 1982; Kenneth, 1986) and of eigenvalue perturbation analysis. The classical framework does not provide the following: (i) the identification of the stable component of Ω \ Γ without root counting at a test point of every component; (ii) the closed-form description of that component as the epigraph of the supremum envelope (Equation 11); and (iii) the global geometric fact that Γ admits no simple counterclockwise loop (Theorem 27), which allows the local right-side instability to be propagated to every non-stable component and thereby yields necessity (see the discussion opening Section 3.4). Items (i)–(iii) are, to our knowledge, new, and they are what turn the classical, essentially computational D-partition picture into the closed-form necessary and sufficient criterion of Theorem 29.

1.2 Organization

The paper is organized as follows. Section 2 introduces the notation, the time-scaling reduction, and the required classical results. Section 3 presents the main results, culminating in the complete stability region characterization (Theorem 29). Section 4 provides three illustrative examples. Section 5 discusses the results in the context of the existing literature and outlines future directions.

2 Preliminaries

Consider the general system (Equation 2) with delay τ > 0. Using the change in variable t ̃ = t / τ, noting thatand that x ( t − τ ) = x ( t ̃ τ − τ ) = x ( ( t ̃ − 1 ) τ ), we obtain (after multiplying by τ 2 and dropping the tildes) the following equation:

d x d t = 1 τ d x d t ̃ , d 2 x d t 2 = 1 τ 2 d 2 x d t ̃ 2 ,

x ̈ t + τ a ̄ 10 x ̇ t + τ 2 a ̄ 00 x t + τ a ̄ 11 x ̇ t − 1 + τ 2 a ̄ 01 x t − 1 = 0 .

Setting a i j = τ 2 − i a ̄ i j, we arrive at the normalized system (Equation 3) with unit delay. From this point on, all stability statements are made for the normalized unit-delay system (Equation 3). The recovery of delay-margin and gain–delay results for general τ is addressed in Remark 11 after the main theorem.

Remark 1 For retarded functional differential equations of the form (Equation 3), asymptotic stability of the trivial solution is equivalent to exponential stability, and both hold if and only if all roots of the characteristic quasi-polynomial (Equation 4) have strictly negative real parts [see, e.g., (Hale and Verduyn Lunel, 1993; Fridman, 2014)]. Moreover, since the solutions of Equation 3 are a time-scaled version of those of Equation 2, the stability of the trivial solution is equivalent for both systems.

The characteristic function of Equation 3 is the quasi-polynomial

H s = h s , e s = e s s 2 + a 10 s + a 00 + a 11 s + a 01 . (4)

It has a principal term s 2 e s [see, e.g. (Ernest Bellman and Kenneth, 1963), for a formal definition]. Writing H ( i w ) = F ( w ) + i G ( w ), the real and imaginary parts are presented in Equation 5:

F w = − a 10 w ⁡ sin w + a 00 − w 2 cos w + a 01 , G w = a 00 − w 2 sin w + a 10 w ⁡ cos w + a 11 w . (5)

2.1 Pontryagin’s theorem

We recall the following classical result, whose necessary condition provides an auxiliary instability certificate used in Section 3.3.

Theorem 1 ([30, Theorem VI]). Let h ( s , e s ) be a quasi-polynomial with a principal term, and let H ( i w ) = F ( w ) + i G ( w ), for real functions F ( w ) and G ( w ). If all zeros of the function H ( s ) have negative real parts, then all zeros of F ( w ) and G ( w ) are real and alternate, as presented in Equation 6:for each w. Furthermore, each of the following conditions is sufficient for all zeros of H ( s ) to have negative real parts:

Remark 2 In the sequel, we only use the necessary implication in Theorem 1, namely, asymptotic stability ⇒ D ( w ) > 0 for all w ∈ R. Although the sufficient conditions (i)–(iii) are classical, they require additional hypotheses on the zeros of F and G that are not easy to verify for our class (Equation 3). Instead, we establish sufficiency through a continuation argument inside the region D (Section 3), which relies on the absence of imaginary-axis roots, path-connectedness of D, and the existence of an explicitly stable anchor point.

Definition 1 (D-partition curve). The auxiliary functions are defined in Equation 7:and the parametric curve Γ : w ↦ A 00 ( w ) , A 10 ( w ), w > 0. By Lemma 2, Γ is the exact locus of imaginary-axis crossings of Equation 4 in the ( a 00 , a 10 )-plane. Therefore, it defines a D-partition curve. The functions A 00, A 10, and the curve Γ will be central to all subsequent developments.

Remark 3 For the second-order case (Equation 3), the function D ( w ) given in Equation 6 admits the factorization shown in Equation 8:which is convenient in proofs since the constraint a 00 = A 00 ( w ) eliminates the second product entirely.

Lemma 2 (Imaginary-axis parametrization). Let A 00 ( w ) and A 10 ( w ) be defined by Equation 7 and H ( s ) be given by Equation 4. For w ≠ 0, That is, the parametric curve Γ : w ↦ A 00 ( w ) , A 10 ( w ) is the exact locus of imaginary-axis crossings in the ( a 00 , a 10 )-plane.

Proof. Writing H ( i w ) = F ( w ) + i G ( w ) = 0 together with Equation 5 amounts to the linear system in the unknowns ( a 00 − w 2 , a 10 ):

cos ⁡ w − w ⁡ sin ⁡ w sin ⁡ w w ⁡ cos ⁡ w a 00 − w 2 a 10 = − a 01 − a 11 w .

The determinant of the first factor is w ( cos 2 ⁡ w + sin 2 ⁡ w ) = w ≠ 0, so the unique solution is a 00 = A 00 ( w ); a 10 = A 10 ( w ).

2.2 Necessary conditions

Theorem 3 ([10, Theorem 3.1]). If the zero solution of Equation 3 is asymptotically stable, Equation 9 can be obtained:

These inequalities define the ambient admissible region Ω (Definition 2 below); outside Ω, stability is impossible.

2.3 Stability regions

Definition 2 Given a pair ( a 01 , a 11 ), let Ω be the region where the necessary conditions (Equation 9) hold Equation 10: For a 00 * ∈ Im A 00, let A 00 − 1 ( a 00 * ) = { w > 0 : A 00 ( w ) = a 00 * } denote the fiber of A 00 at a 00 *. The supremum envelope of the D-partition curve Γ is shown in Equation 11:the supremum is in fact a maximum (Lemma 4). The stability candidate region is presented in Equation 12:

The construction of the supremum envelope is illustrated in Figure 1: for the abscissa a 00 * shown there, the fiber contains three frequencies, w 1, w 2, w 3, and the envelope selects the largest of the corresponding ordinates A 10 ( w i ). Regions Ω and D are illustrated in Figure 2. Since A 00 ( w ) → − a 01 as w → 0 + and A 00 ( w ) → + ∞ as w → ∞, the intermediate value theorem yields ( − a 01 , + ∞ ) ⊂ Im A 00, so A ̂ 10 ( a 00 ) is well defined for every a 00 > − a 01, which is the range required by Ω.

Remark 4 (Regularity of the envelope). Although the proof of Theorem 29 does not require it, it is worth recording the regularity of A ̂ 10. By Lemma 4, each fiber is finite; away from the critical values of A 00 (i.e., for a 00 such that A 00 ′ ( w ) ≠ 0 on the fiber), the implicit function theorem yields finitely many locally real-analytic branches w j ( a 00 ), so A ̂ 10 is piecewise real-analytic, with possible corner points (continuous but non-differentiable) where the maximizing branch switches—these are the junctions between consecutive “scallops” shown in Figure 2b. At a critical value of A 00—a fold of the curve Γ, where A 00 ′ ( w ) = 0—a pair of preimage frequencies is created or annihilated, and A ̂ 10 may exhibit a downward jump discontinuity if the disappearing branch was the maximizing branch. The envelope is, however, always upper semicontinuous: this follows from the closedness of the set S established in the proof of Lemma 11, and it is precisely the property that makes D open. Finally, Corollary 19b shows that, except in the degenerate case β ( w 0 ) = A ̂ 10 ( a 00 ) discussed in Remark 9, the maximum in Equation 11 is attained on an ascending branch of A 00 ( A 00 ′ ( w 0 ) > 0 ).

3 Main results

This section states the main results and summarizes the proof strategy. For fixed ( a 01 , a 11 ), the analysis proceeds in two stages.

Stage 1 (sufficiency by continuation). The imaginary-axis parametrization (Lemma 2) shows that H ( i w ) = 0 occurs precisely on the D-partition curve Γ : w ↦ ( A 00 ( w ) , A 10 ( w ) ). Since every point of D lies strictly above the upper envelope of Γ, no characteristic root can lie on the imaginary axis inside D. Combining this with path-connectedness of D and an explicitly stable anchor point, root continuity implies that every point in D is asymptotically stable.

Stage 2 (necessity by D-partition and topology). Conversely, we show that every point outside D is unstable. Points outside the ambient region Ω violate the classical necessary conditions (Theorem 3). Points on the D-partition curve Γ itself have a purely imaginary characteristic root. For points in Ω \ ( Γ ∪ D ), the argument is topological: the sensitivity identity linking H ′ to Γ ′ (Lemma 6) establishes that points on the right-hand side of Γ are locally unstable (Corollary 8), and the D-partition framework (Corollary 14) propagates this instability to the entire connected component. The key geometric step shows that every such component has at least one right-facing boundary arc, which follows from the exclusion of counterclockwise loops in Γ (Theorem 27); this step is necessary because the classical D-partition framework alone cannot exclude “stable islands” enclosed counterclockwise by Γ (see the discussion opening Section 3.4).

In addition, the crossing-value identity (Section 3.3) provides a complementary Pontryagin-strip instability certificate below the envelope, giving further geometric insight into the boundary structure.

3.1 Imaginary-axis analysis and crossing direction

Lemma 4 For any a 01 , a 11 ∈ R and any a 00 * ∈ R, the fiber A 00 − 1 ( a 00 * ) is finite. In particular, the supremum in Equation 11 is a maximum.

Proof. From Equation 7, A 00 ( w ) = w 2 − a 11 w ⁡ sin ⁡ w − a 01 ⁡ cos ⁡ w. Since | w ⁡ sin ⁡ w | ≤ w and | cos ⁡ w | ≤ 1,so for every fixed a 00 *, there exists W > 0 such that A 00 ( w ) > a 00 * for all w ≥ W. Thus, every solution of A 00 ( w ) = a 00 * lies in ( 0 , W ]. Now, A 00 is real-analytic, being a finite combination of w 2, w ⁡ sin ⁡ w, and cos ⁡ w. If the solution set was infinite, it would accumulate in [ 0 , W ], and the identity theorem would force A 00 ≡ a 00 * on R +, contradicting the unboundedness above.

A 00 w ≥ w 2 − | a 11 | w − | a 01 | → + ∞ w → ∞ ,

Corollary 5 (No imaginary-axis roots inside D). If ( a 00 , a 10 ) ∈ D, then H ( i w ) ≠ 0 for every w ∈ R.

Proof. For w = 0, ( a 00 , a 10 ) ∈ D ⊂ Ω gives H ( 0 ) = a 00 + a 01 > 0. For w ≠ 0, H ( i w ) = 0 would require a 10 = A 10 ( w ) with A 00 ( w ) = a 00 (Lemma 2), hence a 10 ≤ A ̂ 10 ( a 00 ), contradicting the definition of D.

Lemma 6 (Sensitivity identity). Along the D-partition curve Γ, Equation 13 can be obtained: In particular, H ′ ( i w 0 ) ≠ 0 if and only if Γ ′ ( w 0 ) ≠ 0, i.e., every regular point of Γ is a simple imaginary-axis crossing.

Proof. Along Γ, we have H i w , Γ ( w ) = 0 for every w > 0. Writing H ( s , a 00 , a 10 ) = e s ( s 2 + a 10 s + a 00 ) + a 11 s + a 01 to make the parameter dependence explicit, differentiating with respect to w yields

0 = d d w H i w , Γ w = i H ′ i w + ∂ H ∂ a 00 A 00 ′ w + ∂ H ∂ a 10 A 10 ′ w .

Since ∂ H / ∂ a 00 = e s and ∂ H / ∂ a 10 = s e s, evaluating at s = i w yields 0 = i H ′ ( i w ) + e i w A 00 ′ ( w ) + i w e i w A 10 ′ ( w ), whence (Equation 13). For the norms,

Since w 0 > 0, these quantities vanish simultaneously if and only if

A 00 ′ w 0 = A 10 ′ w 0 = 0 .

Corollary 7 (Crossing direction). Fix a 01 , a 11 , a 00 ∈ R with a 00 > − a 01. Let w 0 > 0 satisfy A 00 ( w 0 ) = a 00, and set a 10 0 ≔ A 10 ( w 0 ), so that H ( i w 0 ) = 0 at ( a 00 , a 10 0 ) by Lemma 2. If Γ ′ ( w 0 ) ≠ 0 (equivalently H ′ ( i w 0 ) ≠ 0, by Lemma 6), the implicit function theorem furnishes a root branch s ( a 00 , a 10 ) with s ( a 00 , a 10 0 ) = i w 0, and we can obtain Equation 14:where p = ( a 00 , a 10 0 ) and A 00 ′ ( w 0 ) = 2 w 0 + ( a 01 − a 11 ) sin w 0 − a 11 w 0 ⁡ cos w 0. In particular, if A 00 ′ ( w 0 ) > 0, then increasing a 10 pushes the root into the open left half-plane.

Proof. Since ∂ H / ∂ a 00 = e s and ∂ H / ∂ a 10 = s e s, the implicit function theorem yields

∂ s ∂ a 00 p = − e i w 0 H ′ i w 0 , ∂ s ∂ a 10 p = − i w 0 e i w 0 H ′ i w 0 .

Using Equation 13 to cancel e i w 0, we can obtain the following equation:

∂ s ∂ a 00 p = − 1 i A 00 ′ w 0 − w 0 A 10 ′ w 0 , ∂ s ∂ a 10 p = − i w 0 i A 00 ′ w 0 − w 0 A 10 ′ w 0 .

Taking real parts with | i A 00 ′ − w 0 A 10 ′ | 2 = A 00 ′ 2 + w 0 2 A 10 ′ 2 = | H ′ ( i w 0 ) | 2 yields Equation 14.

Corollary 8 (Right-side instability). Let p = Γ ( w 0 ) be a regular point of Γ ( Γ ′ ( w 0 ) ≠ 0 ). Then, points sufficiently close to p on the right-hand side of the oriented curve Γ are unstable.

Proof. By Lemma 6, H ′ ( i w 0 ) ≠ 0, so the crossing at p is simple, and the implicit function theorem provides a root branch s ( a 00 , a 10 ) with s ( p ) = i w 0. The right unit normal to Γ at p is n R = A 10 ′ ( w 0 ) , − A 00 ′ ( w 0 ) / ‖ Γ ′ ( w 0 ) ‖. By Equation 14,

∇ a 00 , a 10 Re s p ⋅ n R = w 0 A 00 ′ w 0 2 + A 10 ′ w 0 2 | H ′ i w 0 | 2 ‖ Γ ′ w 0 ‖ = w 0 ‖ Γ ′ w 0 ‖ | H ′ i w 0 | 2 > 0 .

Hence, for all sufficiently small ε > 0, the point p + ε n R has a characteristic root with positive real part.

Remark 5 (Degenerate crossings). Lemma 6 shows that the degenerate situation H ′ ( i w 0 ) = 0 occurs precisely at the singular points of the parametrization, Γ ′ ( w 0 ) = 0; there, the imaginary root i w 0 has multiplicity at least two, the crossing may be tangential, and geometrically, Γ exhibits a cusp or a higher-order flat point where crossing frequencies merge. Since A 00 and A 10 are real-analytic and non-constant, such singular points are isolated, so regular points are dense along Γ. Three observations show that these degeneracies do not affect the main result. First, every point of Γ carries a root on the imaginary axis and is therefore not asymptotically stable regardless of the multiplicity or the crossing direction (Case 1 in the proof of Theorem 29). Second, the necessity argument invokes transversality only at one regular boundary point of each component—Lemma 28 uses the density of regular points—and never at all of them. Third, the simultaneous-critical-point configuration A 00 ′ ( w 0 ) = A 10 ′ ( w 0 ) = 0 is analyzed exactly in Lemma 25, and the degenerate envelope maximizer with A 10 ( w 0 ) = β ( w 0 ) is discussed in Remark 9.

Remark 6 Pontryagin’s necessary condition (Theorem 1) requires D ( w ) > 0 for all w ∈ R. The factored form (Equation 8) shows that, at any w 0 with A 00 ( w 0 ) = a 00, D ( w 0 ) vanishes when a 10 = A 10 ( w 0 ). Hence, This yields an alternative (indirect) proof that stability implies membership in D. The present approach establishes this directly through the imaginary-axis parametrization (Lemma 2) and sensitivity analysis (Lemma 6), bypassing D ( w ) entirely.

3.2 Continuation infrastructure

The following four lemmas provide the ingredients for the continuation argument in the sufficiency direction of Theorem 29. They are independent of the Pontryagin factorization and counterclockwise-loop analysis developed in Sections 3.3, Section 3.4 and are collected here so that the final assembly (Section 3.5) can proceed without interruption.

The first lemma provides a uniform bound on the D-partition curve.

Lemma 9 Let A 10 be given by Equation 7 and set M ≔ | a 01 | + | a 11 |. Then, | A 10 ( w ) | ≤ M for all w > 0, and consequently, A ̂ 10 ( a 00 ) ≤ M for every a 00 ∈ Im A 00.

Proof. The triangle inequality yields | A 10 ( w ) | ≤ | a 01 | sin ⁡ w w + | a 11 | | cos ⁡ w | ≤ | a 01 | + | a 11 | = M.

The next lemma establishes that D is path-connected, by showing that it is an epigraph in the a 10-direction.

Lemma 10 The region D is path-connected. Moreover, with M ≔ | a 01 | + | a 11 |, for any ( a 00 , a 10 ) ∈ D, the vertical segment { ( a 00 , a ̃ 10 ) : a ̃ 10 ∈ [ a 10 , M + 2 ] } lies in D. If a 00 ⋆ > − a 01, then the horizontal segment { ( a ̃ 00 , M + 2 ) : a ̃ 00 ∈ [ a 00 , a 00 ⋆ ] } lies in D as well.

Proof. By Definition 2, D is upward closed in a 10: if ( a 00 , a 10 ) ∈ D and a ̃ 10 ≥ a 10, then a ̃ 10 > A ̂ 10 ( a 00 ) and ( a 00 , a ̃ 10 ) ∈ D. By Lemma 9, A ̂ 10 ≤ M; hence, any point in Ω with a 10 = M + 2 belongs to D.

For the horizontal segment at a 10 = M + 2, Lemma 9 yields A ̂ 10 ( a 00 ) ≤ M < M + 2; hence, every point on this segment that lies in Ω also lies in D. It remains to check Ω. The constraints defining Ω in Equation 10 are linear, hence convex; therefore, if the endpoints satisfy them, then the whole segment does. In particular, a 00 ⋆ > − a 01 ensures a 00 + a 01 > 0 along the segment, and a 00 + a 10 + a 11 > 0 holds because a 10 = M + 2 is fixed and positive.

In particular, D is an epigraph in the a 10-direction: its boundary inside Ω is the graph a 10 = A ̂ 10 ( a 00 ), and D is the region above this curve (see Figure 2).

Lemma 11 (D is a connected component of Ω \ Γ). The set D is a connected component of Ω \ Γ.

Proof. By Lemma 10, D is connected. It remains to show that D is both open and closed in Ω \ Γ.

Define

S ≔ a 00 , a 10 ∈ Ω : ∃ w > 0 , A 00 w = a 00 , A 10 w ≥ a 10 .

Then, within Ω \ Γ,

D = Ω \ Γ \ S .

We claim that S is closed in Ω. Let ( a 00 ( n ) , a 10 ( n ) ) ∈ S converge in Ω to ( a 00 * , a 10 * ). Choose w n > 0 such that

A 00 w n = a 00 n , A 10 w n ≥ a 10 n .

Since A 00 ( w ) → + ∞ as w → ∞, the sequence ( w n ) is bounded. Passing to a subsequence if necessary, we may assume w n → w * ≥ 0. Because ( a 00 * , a 10 * ) ∈ Ω, we have a 00 * > − a 01. Since A 00 ( w ) → − a 01 as w → 0 +, it follows that w * > 0. By continuity,so ( a 00 * , a 10 * ) ∈ S. Thus, S is closed in Ω.

A 00 w * = a 00 * , A 10 w * ≥ a 10 * ,

Therefore, D = ( Ω \ Γ ) \ S is open in Ω \ Γ. It is also closed in Ω \ Γ: the boundary of D within Ω consists of points ( a 00 , a 10 ) with a 10 = A ̂ 10 ( a 00 ). By Lemma 4, this maximum is attained at some w * > 0 with A 00 ( w * ) = a 00 and A 10 ( w * ) = a 10, so ( a 00 , a 10 ) = Γ ( w * ) ∈ Γ. Hence, ∂ D ∩ Ω ⊂ Γ, which yields ∂ D ∩ ( Ω \ Γ ) = ∅. Therefore, D is closed in Ω \ Γ.

Thus, D is clopen and connected in Ω \ Γ, so it is a connected component of Ω \ Γ.

The following lemma provides an anchor point where stability can be verified by elementary means (Rouché’s theorem).

Lemma 12 Let M ≔ | a 01 | + | a 11 |. If a 10 > M and a 00 > 1 + M, then all zeros of H ( s ) satisfy R e ( s ) < 0.

Proof. Since e s ≠ 0, the zeros of H coincide with those of Δ ( s ) ≔ e − s ⁡ H ( s ) = P ( s ) + Q ( s ), where P ( s ) = s 2 + a 10 s + a 00 and Q ( s ) = ( a 11 s + a 01 ) e − s. Since a 10 > 0 and a 00 > 0, P is Hurwitz. We verify | Q | < | P | on the boundary of Σ R ≔ { s : R e ( s ) > 0 , | s | < R } for all sufficiently large R. On s = i ω, | Q ( i ω ) | ≤ | a 11 ‖ ω | + | a 01 |. For | ω | ≥ 1, this yields | Q | ≤ M | ω | < a 10 | ω | ≤ | P |. For | ω | ≤ 1, | Q | ≤ M < a 00 − 1 ≤ a 00 − ω 2 ≤ | P |. On the arc | s | = R, R e ( s ) ≥ 0: | Q | ≤ M R, while | P | ≥ R 2 − a 10 R − a 00 > M R when R is sufficiently large. By Rouché’s theorem, Δ and P have the same number of zeros in Σ R, namely, zero, for all sufficiently large R. Thus, Δ has no zeros in the closed right half-plane, and therefore, all zeros of H satisfy R e ( s ) < 0.

The last lemma of this subsection is the root-continuity principle. It will also be used directly in the necessity argument (via Corollary 14).

Lemma 13 Along any continuous path γ : [ 0,1 ] → R 4 in parameter space ( a 00 , a 10 , a 01 , a 11 ) such that H ( i ω ) ≠ 0 for all ω ∈ R and all parameter values on γ, the number of zeros in Equation 4 in the open right half-plane (counted with multiplicity) is constant along γ.

Proof. Set K ≔ γ ( [ 0,1 ] ) and choose B > 0 such that | a i j | ≤ B on K. For λ = ( a 00 , a 10 , a 01 , a 11 ) ∈ K, write Δ λ ( s ) ≔ e − s H λ ( s ) = s 2 + a 10 s + a 00 + ( a 11 s + a 01 ) e − s .. Since e − s ≠ 0, H λ and Δ λ have the same zeros. For R e s ≥ 0,

Hence, there exists R > 0, independent of λ ∈ K, such that Δ λ ( s ) ≠ 0 whenever R e s ≥ 0 and | s | ≥ R.

Therefore, all zeros of H λ in the open right half-plane lie in

Σ R ≔ s : R e s > 0 , | s | < R .

On the arc { | s | = R , R e s ≥ 0 }, there are no zeros by construction, and on the segment { i ω : | ω | ≤ R }, there are none by hypothesis. Thus, H λ has no zeros on ∂ Σ R for any λ ∈ K. By the argument principle, the numberequals the number of zeros of H λ in Σ R, counted with multiplicity. Since the integrand depends continuously on λ and N ( λ ) ∈ Z, the function N is constant on the connected set K. Hence, the number of zeros in Equation 4 in the open right half-plane remains constant along γ.

N λ = 1 2 π i ∫ ∂ Σ R H λ ′ s H λ s d s

Corollary 14 (D-partition property). For fixed ( a 01 , a 11 ), the number of zeros in Equation 4 in the open right half-plane is constant on each connected component of Ω \ Γ.

Proof. Let γ : [ 0,1 ] → Ω \ Γ be a continuous path. For every ( a 00 , a 10 ) = γ ( t ), H ( 0 ) = a 00 + a 01 > 0 because γ ( t ) ∈ Ω, and H ( i w ) ≠ 0 for w ≠ 0 by Lemma 2 because γ ( t ) ∉ Γ. By Lemma 13, the number of zeros in the open right half-plane is constant along γ, hence constant on each path-connected component of Ω \ Γ. Since Ω is open and Γ ∩ Ω is relatively closed in Ω, the set Ω \ Γ is open, hence locally path-connected, and its connected and path-connected components coincide.

Having established all the ingredients, we can now package sufficiency as a self-contained result.

Proposition 15 (Sufficiency of D). Every ( a 00 , a 10 ) ∈ D is asymptotically stable.

Proof. Set M ≔ | a 01 | + | a 11 | and define

q ≔ M + 2 , M + 2 .

Since M ≥ | a 01 | and M ≥ | a 11 |,so q ∈ Ω. Moreover, by Lemma 9,hence q ∈ D. Since also M + 2 > M and M + 2 > 1 + M, Lemma 12 shows that the system at q is asymptotically stable.

M + 2 + a 01 ≥ 2 > 0 , M + 2 + M + 2 + a 11 ≥ M + 4 > 0 ,

A ̂ 10 M + 2 ≤ M < M + 2 ,

Now, let ( a 00 , a 10 ) ∈ D be arbitrary. By Lemma 10, there exists a continuous path in D joining ( a 00 , a 10 ) to q. Along this path, Corollary 5 yieldsTherefore, by Lemma 13, the number of zeros in Equation 4 in the open right half-plane remains constant along the path. Since this number is zero at q, it is also zero at ( a 00 , a 10 ).

H i ω ≠ 0 for all ω ∈ R .

Finally, Corollary 5 excludes zeros on the imaginary axis at ( a 00 , a 10 ). Hence, all zeros of Equation 4 satisfy R e s < 0, and the system is asymptotically stable.

3.3 Crossing-value identity and Pontryagin factorization

The main necessity argument (Section 3.5) relies on D-partition topology and the right-side instability criterion of Corollary 8. Here, we develop a complementary algebraic route. Lemma 16 yields the identitywhich introduces the auxiliary function β and leads to a factorization of the Pontryagin function at constrained frequencies (Lemma 17). Together with Pontryagin’s necessary condition, this yields an explicit unstable stripimmediately below the envelope (Corollary 19c). The same identity also determines the sign of A 00 ′ ( w 0 ) at an envelope maximizer, whenever β ( w 0 ) < A ̂ 10 ( a 00 ) (Corollary 19b).

A 00 ′ w = w A 10 w − β w ,

β w 0 < a 10 < A ̂ 10 a 00

Lemma 16 (Crossing-value identity). For every w > 0, Equivalently, ifthen Equation 16 can be obtained: In particular, sgn A 00 ′ ( w ) = sgn A 10 ( w ) − β ( w ) for every w > 0.

Proof. Differentiating A 00 ( w ) in Equation 7 yields

A 00 ′ w = 2 w + a 01 − a 11 sin ⁡ w − a 11 w ⁡ cos ⁡ w .

On the other hand, from the expression of A 10 ( w ) in Equation 7, we obtain

w 2 + A 10 w = 2 w + a 01 ⁡ sin ⁡ w − a 11 w ⁡ cos ⁡ w .

Subtracting the second identity from the first yields Equation 15. Since w β ( w ) = a 11 ⁡ sin ⁡ w − 2 w, Equation 16 follows by substitution.

Lemma 17 (Pontryagin factorization at constrained frequencies). Let w 0 > 0 be satisfied A 00 ( w 0 ) = a 00, and set v 0 ≔ A 10 ( w 0 ) and β 0 ≔ β ( w 0 ) = a 11 ⁡ sin ( w 0 ) / w 0 − 2. Then, we can obtain Equation 17: Since w 0 2 > 0, D ( w 0 ) < 0 if and only if a 10 lies in the open interval ( min ( v 0 , β 0 ) , max ( v 0 , β 0 ) ).

Proof. Since A 00 ( w 0 ) = a 00, the second product in Equation 8 vanishes, yielding

D w 0 = v 0 − a 10 a 11 w 0 ⁡ sin w 0 − w 0 2 a 10 + 2 .

The second factor rewrites as w 0 2 a 11 ⁡ sin w 0 w 0 − a 10 − 2 = w 0 2 ( β 0 − a 10 ), yielding Equation 17. The sign criterion follows from the factorization and w 0 2 > 0.

The following corollary extends Lemma 17 to arbitrary frequencies. Although it is not used directly in the main proof, the decomposition provides useful geometric insight: for each fixed w > 0, the inequality D ( w ) < 0 can occur only if a 10 lies between A 10 ( w ) and β ( w ), and even then, only when a 00 lies in a bounded strip around A 00 ( w ).

Corollary 18 (Signed decomposition of D ( w )). For every w > 0 and every ( a 00 , a 10 ) ∈ R 2, we can obtain Equation 18:

Proof. The two brackets in the factored form (Equation 8) simplify as follows. For the first, a 11 w ⁡ sin ⁡ w − w 2 ( a 10 + 2 ) = w 2 a 11 ⁡ sin ⁡ w w − 2 − a 10 = w 2 ( β ( w ) − a 10 ). For the second, w 2 − a 00 − a 01 ⁡ cos ⁡ w − a 11 w ⁡ sin ⁡ w = ( w 2 − a 11 w ⁡ sin ⁡ w − a 01 ⁡ cos ⁡ w ) − a 00 = A 00 ( w ) − a 00. Substituting into Equation 8 yields Equation 18.

Remark 7 The decomposition in Equation 18 separates D ( w ) into a signed term and a nonnegative term: The square term vanishes precisely at the preimage frequencies A 00 ( w ) = a 00, recovering Lemma 17. The first term is a signed quadratic in a 10, with roots A 10 ( w ) and β ( w ); in particular, it is negative if and only if a 10 lies strictly between A 10 ( w ) and β ( w ). Hence, D ( w ) < 0 requires

Corollary 19 (Ascending branch at the envelope; Pontryagin strip). Fix a 01 , a 11 ∈ R and a 00 > − a 01 , and let w 0 > 0 be a maximizer with A 00 ( w 0 ) = a 00 and v 0 ≔ A 10 ( w 0 ) = A ̂ 10 ( a 00 ) . Set β 0 ≔ β ( w 0 ). By Theorem 1, the system is not asymptotically stable.

Proof. (a) For a 10 > A ̂ 10 ( a 00 ), ( a 00 , a 10 ) ∈ D, hence asymptotically stable by Proposition 15, so D ( w 0 ) > 0 by Theorem 1. In Equation 17, the first factor v 0 − a 10 = A ̂ 10 ( a 00 ) − a 10 < 0 and w 0 2 > 0, hence necessarily β 0 − a 10 < 0. Thus, β 0 < a 10 for every a 10 > A ̂ 10 ( a 00 ), and therefore, β 0 ≤ A ̂ 10 ( a 00 ).

(b) By Equation 16,since w 0 > 0 and v 0 − β 0 = A ̂ 10 ( a 00 ) − β 0 > 0.

A 00 ′ w 0 = w 0 v 0 − β 0 > 0 ,

(c) If a 10 ∈ ( β 0 , A ̂ 10 ( a 00 ) ), then v 0 − a 10 > 0, β 0 − a 10 < 0, and w 0 2 > 0. Hence, Equation 17 yields D ( w 0 ) < 0, contradicting the necessary condition for asymptotic stability in Theorem 1.

Remark 8 Corollary 19b shows that A 00 ′ ( w 0 ) > 0 at an envelope maximizer need not be assumed separately: it follows from the identity Equation 16 whenever β ( w 0 ) < A ̂ 10 ( a 00 ).

Lemma 20 (Boundary nonnegativity). Fix a 00 > − a 01 and set v ≔ A ̂ 10 ( a 00 ). Then,where D ( w ) is evaluated at ( a 00 , a 10 ) = ( a 00 , v ).

Proof. For every ε > 0, ( a 00 , v + ε ) ∈ D, hence asymptotically stable by Proposition 15. Pontryagin’s necessary condition yields D ( w ) > 0 for all w > 0 at ( a 00 , v + ε ). Since D ( w ) depends continuously on a 10, letting ε ↓ 0 yields D ( w ) ≥ 0 at ( a 00 , v ).

Remark 9 If w 0 is an envelope maximizer with A 10 ( w 0 ) = β ( w 0 ) = A ̂ 10 ( a 00 ), then A 00 ′ ( w 0 ) = 0 by Equation 16. Since A 00 ( w 0 ) = a 00 and A 00 ′ ( w 0 ) = 0, the square term in Equation 18 is O ( w − w 0 ) 4, and the local expansion of D near w 0 yields Lemma 20 then forces A 10 ′ ( w 0 ) β ′ ( w 0 ) ≥ 0. This rules out part of the degenerate case but does not exclude it completely.

3.4 Exclusion of counterclockwise simple loops

Before entering the technical development, let us explain why this subsection is needed at all—that is, why the necessity part of Theorem 29 cannot be concluded directly from the ingredients already available. The D-partition property (Corollary 14) shows that the number of unstable roots is constant on every connected component of Ω \ Γ, and the right-side instability criterion (Corollary 8) produces unstable points next to every regular arc of Γ—but only on its right-hand side. The missing step is topological: a priori, some component C ≠ D could be bounded exclusively by arcs having C on their left-hand side, i.e., Γ could enclose C by a simple loop traversed counterclockwise. For such a hypothetical “stable island” below the envelope, Corollary 8 would place unstable points outside C, and no instability certificate would propagate into it. This configuration cannot be excluded by the D-partition machinery alone—it occurs for general parametric curves—and, if it were to occur here, the region D would fail to be the complete stability region. The purpose of this subsection is to prove that the specific curve Γ generated by Equation 7 admits no simple counterclockwise loop (Theorem 27); Lemma 28 then guarantees that every component other than D has at least one right-facing boundary arc, through which the instability enters. This is the only step of the proof that is not part of the classical D-partition toolkit, and it is where the specific analytic structure of A 00 and A 10 enters.

We prove that the parametric curvecannot generate a simple counterclockwise loop. The argument proceeds in two steps: Proposition 22 excludes loops whose highest point is an interior regular point, and Proposition 26 excludes “teardrop” loops whose highest point is the closing point. We begin by establishing the algebraic formulae needed for both cases.

Γ w = A 00 w , A 10 w , w > 0 ,

Lemma 21 (Critical points of A 10 away from tan ⁡ w = w). Let w > 0 be satisfied A 10 ′ ( w ) = 0 and δ ( w ) ≔ w ⁡ cos ⁡ w − sin ⁡ w ≠ 0 . Set κ ( w ) ≔ w 2 − sin 2 ⁡ w and μ ( w ) ≔ w 2 + w ⁡ sin ⁡ w ⁡ cos ⁡ w − 2 ⁡ sin 2 ⁡ w . Then, a. b. c. Moreover, κ ( w ) > 0 and μ ( w ) > 0 for all w > 0.

Proof. Setting A 10 ′ ( w ) = 0 inyields (a).

A 10 ′ w = a 01 w ⁡ cos ⁡ w − sin ⁡ w w 2 + a 11 ⁡ sin ⁡ w

Substituting this intoyieldswhich is (b).

A 00 ′ w = 2 w + a 01 − a 11 sin ⁡ w − a 11 w ⁡ cos ⁡ w

A 00 ′ w = 2 w − a 11 w 2 ⁡ sin 2 ⁡ w δ w − a 11 ⁡ sin ⁡ w − a 11 w ⁡ cos ⁡ w = 2 w δ w − a 11 w 2 ⁡ sin 2 ⁡ w + sin ⁡ w δ w + w ⁡ cos ⁡ w δ w δ w = − a 11 κ w + 2 w 2 ⁡ cos ⁡ w − w ⁡ sin ⁡ w δ w ,

Similarly, substituting (a) intoyields (c) after simplification.

A 10 ″ w = a 01 − sin ⁡ w w − 2 ⁡ cos ⁡ w w 2 + 2 ⁡ sin ⁡ w w 3 + a 11 ⁡ cos ⁡ w

For the positivity claims,since | sin ⁡ w | < w for w > 0.

κ w = w 2 − sin 2 ⁡ w > 0 w > 0 ,

For μ, note that μ ( 0 ) = 0, μ ′ ( 0 ) = 0, and

μ ′ ′ w = 4 ⁡ sin ⁡ w sin ⁡ w − w ⁡ cos ⁡ w .

On ( 0 , π ) both factors are positive: sin ⁡ w > 0, andwhile sin ⁡ 0 − 0 ⁡ cos ⁡ 0 = 0. Hence, μ ′ ′ ( w ) > 0 on ( 0 , π ). Since μ ′ ( 0 ) = 0, it follows that μ ′ ( w ) > 0 on ( 0 , π ); since μ ( 0 ) = 0, we conclude that μ ( w ) > 0 on ( 0 , π ). For w ≥ 2,

d d w sin ⁡ w − w ⁡ cos ⁡ w = w ⁡ sin ⁡ w > 0 ,

μ w ≥ w 2 − w 2 − 2 > 0 .

Therefore, μ ( w ) > 0 for all w > 0.

Remark 10 The case δ ( w ) = 0 (i.e., tan ⁡ w = w) with A 10 ′ ( w ) = 0 forces a 11 = 0, since thenand sin ⁡ w ≠ 0 whenever tan ⁡ w = w with w > 0. The formulae in Lemma 21 do not apply in this degenerate case.

Proposition 22 (No counterclockwise loop with the regular top point). The curve Γ admits no simple counterclockwise loop whose highest point is attained at an interior regular point of the parametrization.

Proof. Suppose, for contradiction, that Γ | [ u , v ] is a simple counterclockwise loop and that w ∗ ∈ ( u , v ) is an interior regular point where A 10 attains its maximum along the loop. Then,and the tangent at the top of a counterclockwise loop points left A 00 ′ ( w ∗ ) < 0.

A 10 ′ w ∗ = 0 , A 10 ′ ′ w ∗ ≤ 0 ,

Case δ ( w ∗ ) ≠ 0. Lemma 21 applies. Since 2 ( w ∗ 2 ⁡ cos w ∗ − w ∗ ⁡ sin w ∗ ) = 2 w ∗ δ ( w ∗ ), formula (b) yields A 00 ′ ( w ∗ ) = ( − a 11 κ + 2 w ∗ δ ) / δ, where κ = w ∗ 2 − sin 2 w ∗ > 0. The inequality A 00 ′ ( w ∗ ) < 0 then rearranges to a 11 / δ > 2 w ∗ / κ > 0. By formula (c), A 10 ″ ( w ∗ ) = a 11 μ ( w ∗ ) / ( w ∗ δ ), and since μ ( w ∗ ) > 0, w ∗ > 0, and a 11 / δ > 0, we obtain A 10 ″ ( w ∗ ) > 0, contradicting A 10 ″ ( w ∗ ) ≤ 0.

Case δ ( w ∗ ) = 0. As noted in the remark following Lemma 21, a 11 = 0. Then, A 10 ( w ) = a 01 ⁡ sin ⁡ w / w, so A 10 ′ ′ ( w ∗ ) = − a 01 ⁡ sin w ∗ / w ∗ and A 00 ′ ( w ∗ ) = 2 w ∗ + a 01 ⁡ sin w ∗. The condition A 10 ′ ′ ( w ∗ ) ≤ 0 requires a 01 ⁡ sin w ∗ ≥ 0, while A 00 ′ ( w ∗ ) < 0 requires a 01 ⁡ sin w ∗ < − 2 w ∗ < 0—a contradiction.

Lemma 23 (Vertical tangencies with the nonnegative second derivative). Let w > 0 be satisfied DefineThen, a. b. If additionally A 00 ″ ( w ) ≥ 0, then

Proof. The positivity of K follows from | sin ⁡ w | < w for w > 0. In addition,so in fact

L 0 = 0 , L ′ w = 2 ⁡ sin 2 ⁡ w ≥ 0 ,

L w = ∫ 0 w 2 ⁡ sin 2 ⁡ t d t > 0 w > 0 .

Assume now that A 00 ′ ( w ) = 0.

First consider the case sin ⁡ w ≠ 0. Solvingfor a 01 yields

2 w + a 01 − a 11 sin ⁡ w − a 11 w ⁡ cos ⁡ w = 0

a 01 = a 11 1 + w ⁡ cot ⁡ w − 2 w ⁡ csc ⁡ w .

Substituting into A 10 ′ and A 00 ″ yieldswhere

A 10 ′ w = a 11 K − 2 w δ w 2 ⁡ sin ⁡ w , A 00 ″ w = a 11 L − 2 δ sin ⁡ w ,

δ ≔ w ⁡ cos ⁡ w − sin ⁡ w .

We also note the identity

K = w L + sin ⁡ w δ .

If sin ⁡ w > 0, then A 00 ″ ( w ) ≥ 0 implies

a 11 ≥ 2 δ L .

Hence,

a 11 K − 2 w δ ≥ 2 δ L K − 2 w δ = 2 δ K − w L L = 2 ⁡ sin ⁡ w δ 2 L ≥ 0 .

Since w 2 ⁡ sin ⁡ w > 0, it follows that A 10 ′ ( w ) ≥ 0.

If sin ⁡ w < 0, then A 00 ″ ( w ) ≥ 0 implies

a 11 ≤ 2 δ L .

Thus,

a 11 K − 2 w δ ≤ 2 δ L K − 2 w δ = 2 ⁡ sin ⁡ w δ 2 L ≤ 0 .

Since now w 2 ⁡ sin ⁡ w < 0, we again conclude that A 10 ′ ( w ) ≥ 0.

Finally, consider the case sin ⁡ w = 0, i.e., w = n π for some positive integer n. Then, cos ⁡ w ≕ c = ± 1, and the condition A 00 ′ ( w ) = 0 reduces toso

2 w − a 11 w c = 0 ,

a 11 c = 2 .

Using the original formulawe obtain at sin ⁡ w = 0

A 00 ″ w = 2 + a 01 ⁡ cos ⁡ w − 2 a 11 ⁡ cos ⁡ w + a 11 w ⁡ sin ⁡ w ,

A 00 ″ w = 2 + a 01 c − 2 a 11 c = a 01 c − 2 .

Hence, A 00 ″ ( w ) ≥ 0 implies

a 01 c ≥ 2 .

However,

A 10 ′ w = a 01 w ⁡ cos ⁡ w − sin ⁡ w w 2 + a 11 ⁡ sin ⁡ w = a 01 c w ≥ 2 w > 0 .

Therefore, in all cases,

A 00 ′ w = 0 , A 00 ″ w ≥ 0 ⇒ A 10 ′ w ≥ 0 .

Lemma 24 If a 11 = 0, the curve Γ admits no simple counterclockwise loop.

Proof. When a 11 = 0, the identity in Equation 16 yields A 00 ′ ( w ) = w A 10 ( w ) + 2 for all w > 0. Suppose, for contradiction, that Γ | [ u , v ] is a simple counterclockwise loop, and set h ≔ max [ u , v ] A 10.

Case h ≤ − 2. Then, A 10 ( w ) + 2 ≤ 0 on [ u , v ], so A 00 ′ ( w ) = w ( A 10 ( w ) + 2 ) ≤ 0, with equality only on the finite zero set of the real-analytic function A 10 + 2. Hence, A 00 is strictly decreasing on [ u , v ], contradicting A 00 ( u ) = A 00 ( v ).

Case h > − 2 , maximum achieved in the interior. Let w ∗ ∈ ( u , v ) be achieved the maximum h. For c ∈ ( − 2 , h ) sufficiently close to h, the connected component of { w ∈ ( u , v ) : A 10 ( w ) > c } containing w ∗ is a single interval ( w − , w + ) ⊂ ( u , v ). On this interval, A 10 ( w ) > c > − 2, so A 00 ′ ( w ) > 0, and A 00 is strictly increasing: the top arc is traversed from left to right. On any simple counterclockwise loop, however, the topmost arc is traversed from right to left, representing a contradiction.

Case h > − 2 , maximum achieved only at the endpoints. Then, A 10 ( u ) = A 10 ( v ) = h and A 10 ( w ) < h on ( u , v ). Since h > − 2 and u , v > 0,

A 00 ′ u = u h + 2 > 0 , A 00 ′ v = v h + 2 > 0 .

Hence, for w > u sufficiently close to u, the branch moves to the right, while A 10 ( w ) < h shows that it lies strictly below the closing point. Thus, the loop leaves the closing point into the lower-right neighborhood. Similarly, for w < v sufficiently close to v, the branch is again to the right of the closing point and strictly below it, so it approaches the closing point from the lower-left neighborhood. This local orientation is clockwise, contradicting the assumption that the loop is counterclockwise.

Lemma 25 (Sign of A 00 ″ at simultaneous critical points). Let a 11 ≠ 0 and suppose that w > 0 satisfies A 00 ′ ( w ) = 0 and A 10 ′ ( w ) = 0. Then, A 00 ″ ( w ) < 0.

Proof. Write s = sin ⁡ w, c = cos ⁡ w, δ = w c − s, and K = w 2 − s 2. By Remark 10, δ ≠ 0, so Lemma 21 applies. From A 10 ′ ( w ) = 0: a 01 = − a 11 w 2 ⁡ s / δ. From A 00 ′ ( w ) = ( − a 11 K + 2 w δ ) / δ = 0: a 11 = 2 w δ / K (since K > 0). Substituting into a 01 yields a 01 = − ( 2 w δ / K ) ( w 2 s / δ ) = − 2 w 3 s / K.

Now, compute A 00 ″ ( w ) = 2 + a 01 ⁡ cos ⁡ w − 2 a 11 ⁡ cos ⁡ w + a 11 w ⁡ sin ⁡ w:

A 00 ″ w = 2 + − 2 w 3 s c K − 4 w δ c K + 2 w 2 s δ K = 2 + 2 w − w 2 s c − 2 δ c + w s δ K .

Expanding δ = w c − s in the numerator,

− w 2 s c − 2 δ c + w s δ = − w 2 s c − 2 w c − s c + w s w c − s = − w 2 s c − 2 w c 2 + 2 s c + w 2 s c − w s 2 = 2 s c − w 2 c 2 + s 2 = 2 s c − w 1 + c 2 .

Therefore,

K 2 A 00 ″ w = K + w 2 s c − w 1 + c 2 = w 2 − s 2 + 2 w s c − w 2 − w 2 c 2 .

Collecting terms,

K 2 A 00 ″ w = − s 2 + 2 w s c − w 2 c 2 = − s − w c 2 = − δ 2 .

Since δ ≠ 0, we have δ 2 > 0 and K > 0; hence, A 00 ″ ( w ) = − 2 δ 2 / K < 0.

Proposition 26 (No counterclockwise teardrop loop). The curve Γ ( w ) = ( A 00 ( w ) , A 10 ( w ) ), w > 0, admits no simple counterclockwise loop whose highest point is the closing point itself.

Proof. The argument proceeds in three steps: first, we locate the first stationary point w m of A 10 after leaving the top and show that A 00 ′ ( w m ) > 0; next, we find the first zero w L of A 00 ′ before w m; finally, we use the vertical-tangency lemma at w L to obtain a contradiction.

Since a 11 = 0 is already excluded by Lemma 24, we may assume a 11 ≠ 0.

Assume, for contradiction, that Γ | [ u , v ] is a simple counterclockwise loop whose maximum of A 10 is attained at the closing point

Γ u = Γ v .

Immediately after leaving the closing point along one branch, the loop moves downward and to the left. Hence, near u,

A 10 ′ < 0 , A 00 ′ < 0 .

Since the loop is closed with A 10 ( u ) = A 10 ( v ), the function A 10 decreases from A 10 ( u ) and must later return to the same value. Therefore, A 10 ′ has at least one zero in ( u , v ). Let w m ∈ ( u , v ) be the first such zero:

A 10 ′ w m = 0 , A 10 ′ w < 0 u < w < w m .

We now show A 00 ′ ( w m ) > 0. If A 00 ′ ( w m ) ≠ 0, since w m is a local minimum of A 10 | [ u , v ] and the loop is simple counterclockwise, the tangent direction at w m must point to the right, yielding A 00 ′ ( w m ) > 0.

If instead A 00 ′ ( w m ) = 0, then Γ ′ ( w m ) = 0. By Lemma 25, A 00 ″ ( w m ) < 0. However, A 00 ′ < 0 on ( u , w m ) and A 00 ′ ( w m ) = 0, so A 00 ′ is increasing at w m, yielding A 00 ″ ( w m ) ≥ 0. This contradicts A 00 ″ ( w m ) < 0. Hence, A 00 ′ ( w m ) = 0 is impossible.

Therefore, A 00 ′ ( w m ) > 0. However, near u, we have A 00 ′ < 0, so by continuity, there exists a first pointsuch that

w L ∈ u , w m

A 00 ′ w L = 0 .

Because w L is the first zero of A 00 ′ after leaving a region where A 00 ′ < 0, we necessarily have

A 00 ″ w L ≥ 0 .

Moreover, since w L < w m, we are still in the interval where

A 10 ′ w L < 0 .

This contradicts Lemma 23, which yields

A 00 ′ w L = 0 , A 00 ″ w L ≥ 0 ⇒ A 10 ′ w L ≥ 0 .

Therefore, no such counterclockwise teardrop loop can exist.

Theorem 27 (No simple counterclockwise loops). The parametric curve Γ ( w ) = ( A 00 ( w ) , A 10 ( w ) ), w > 0, admits no simple counterclockwise loop.

Proof. If a 11 = 0, the claim follows from Lemma 24. Assume a 11 ≠ 0. By Proposition 22, a simple counterclockwise loop cannot have an interior regular top point. Hence, the only remaining possibility would be a teardrop loop whose highest point is the closing point itself. This case is excluded by Proposition 26. Therefore, no simple counterclockwise loop exists.

3.5 Necessity and sufficiency of D for stability

We now assemble the previously established ingredients to prove that D is exactly the stability region.

Recall that the D-partition boundary is the parametric curveand that the upper envelope is

Γ w = A 00 w , A 10 w , w > 0 ,

A ̂ 10 a 00 ≔ sup A 10 w : A 00 w = a 00 .

Accordingly,

D ≔ a 00 , a 10 ∈ Ω : a 10 > A ̂ 10 a 00 .

We shall use the following facts established in the previous subsections:

• F1. If ( a 00 , a 10 ) ∈ D, then the system is asymptotically stable (Proposition 15).
• F2. The curve Γ is a D-partition boundary: the number N RHP of characteristic roots in the open right half-plane is constant on each connected component of Ω \ Γ (Corollary 14).
• F3. At every regular point of Γ, points sufficiently close on the right normal ray are unstable (Corollary 8).
• F4. The curve Γ admits no simple counterclockwise loop (Theorem 27).

The only remaining step is to convert F4 into the topological statement needed to propagate the local right-side instability to every component outside D.

The following lemma is the sole topological use of the no-counterclockwise-loop theorem: if every boundary arc of a component C had C on the left, one could extract a simple counterclockwise loop in Γ.

Lemma 28 (Every non-D component has a right-facing boundary arc). Let C be a connected component of Ω \ Γ with C ≠ D. Then, there exists a regular point,such that C lies locally on the right-hand side of the oriented curve Γ at p.

Proof. Assume, for contradiction, that at every regular boundary point p ∈ ∂ C ∩ Γ, the component C lies locally on the left-hand side of the oriented curve.

First, ∂ C ∩ Γ ≠ ∅: since C ≠ D and both C , D are connected components of Ω \ Γ contained in the connected set Ω, any path in Ω from C to D must cross Γ, so ∂ C meets Γ. Moreover, since Γ is real-analytic, its singular points are isolated, and regular points are dense in ∂ C ∩ Γ.

By Lemma 9, | A 10 ( w ) | ≤ M for all w > 0, so Γ is confined to the horizontal strip | a 10 | ≤ M. Choose R > 0 large enough that the disk B R of radius R, centered at the origin, contains every bounded connected component of Γ − 1 ( B R ) that intersects ∂ C, and consider C R ≔ C ∩ B R. Since C R is bounded and open, its boundary within Ω consists of arcs on Γ and arcs on the circle ∂ B R. The boundary arcs on Γ are nonempty for R sufficiently large (since ∂ C ∩ Γ ≠ ∅).

Because C is always locally on the left of Γ at every regular boundary point, the Γ-arcs forming ∂ C R are oriented so that C R lies on their left. These arcs, together with arcs on ∂ B R, form a closed boundary cycle of the bounded region C R. By removing self-intersecting sub-arcs if necessary (using the fact that Γ is real-analytic with isolated self-intersections), one extracts a simple closed sub-curve of Γ that is oriented counterclockwise with respect to the region it encloses. This contradicts Fact F4.

Therefore, at least one regular boundary point of C must have C locally on the right-hand side of Γ.

We can now prove the main result.

Theorem 29 (Stability region). Given a 01 , a 11 ∈ R, the zero solution of Equation 3 is asymptotically stable if and only if ( a 00 , a 10 ) ∈ D, where D is defined by Equation 2.

Proof. (Sufficiency). Immediate from Proposition 15.

(Necessity). Let ( a 00 , a 10 ) ∉ D.

Case 0: ( a 00 , a 10 ) ∉ Ω. Then, at least one of the necessary conditions (Equation 9) is violated, so the system is not asymptotically stable by Theorem 3.

Case 1: ( a 00 , a 10 ) ∈ Γ. Then, there exists w 0 > 0 withso H ( i w 0 ) = 0 by Lemma 2. Hence, the system is not asymptotically stable.

a 00 = A 00 w 0 , a 10 = A 10 w 0 ,

Case 2: ( a 00 , a 10 ) ∈ Ω \ ( Γ ∪ D ). Let C be the connected component of Ω \ Γ containing ( a 00 , a 10 ). Then, C ≠ D since D is a connected component of Ω \ Γ by Lemma 11. By Lemma 28, there exists a regular pointsuch that C lies locally on the right-hand side of Γ at p.

p ∈ ∂ C ∩ Γ

By Fact F3, there exists ε 0 > 0 such that for every sufficiently small ε ∈ ( 0 , ε 0 ), the pointlies in C and is unstable. Hence, N RHP > 0 at one point of C. By Fact F2, N RHP is constant throughout the whole connected component C. Therefore, every point of C is unstable; in particular, ( a 00 , a 10 ) is unstable.

p + ε n R p

Thus, every point outside D is non-stable, whereas every point inside D is stable. Hence,

a 00 , a 10 ∈ D ⟺ the system is asymptotically stable.

Remark 11 (Recovering delay margins and gain–delay charts). Consider the same dynamics with an explicit delay τ > 0, as presented in Equation 19:The characteristic equation of Equation 19 is obtained from Equation 3 by replacing e − s with e − s τ. With the time scaling t = τ r and λ = τ s, one obtains the unit-delay form of Equation 3 with scaled coefficients Therefore, Theorem 29 can be applied for each τ by evaluating the region D associated with ( a 01 , a 11 ) = ( τ 2 a ̄ 01 , τ a ̄ 11 ) at the point ( a 00 , a 10 ) = ( τ 2 a ̄ 00 , τ a ̄ 10 ). In particular, the delay margin can be computed as the largest τ such thatwhere A ̂ 10 ( a 01 , a 11 ) denotes the supremum envelope for the fixed pair ( a 01 , a 11 ). This reduces delay-margin computations to a one-dimensional search over τ, where each evaluation can be carried out with Algorithm 1.

3.6 Algorithmic computation of the stability boundary

For fixed ( a 01 , a 11 ), the boundary of D inside Ω is the graphwhere A ̂ 10 is the supremum envelope (Equation 11). By Lemma 2, the D-partition curve Γ : w ↦ ( A 00 ( w ) , A 10 ( w ) ) is the exact imaginary-axis crossing locus, and the boundary of D is the upper envelope of Γ.

a 10 = A ̂ 10 a 00 , a 00 > − a 01 ,

The search for preimage points w satisfying A 00 ( w ) = a 00 can be confined to a finite interval:

Corollary 30 All solutions of A 00 ( w ) = a 00 in R > 0 lie in the interval [ 0 , W ( a 00 ) ], where

Proof. From A 00 ( w ) ≥ w 2 − | a 11 | w − | a 01 | (as in the proof of Lemma 4), A 00 ( w ) > a 00 holds whenever w 2 − | a 11 | w − | a 01 | > a 00, i.e., whenever w exceeds the positive root of the quadratic w 2 − | a 11 | w − ( | a 00 | + | a 01 | ) = 0. The stated W ( a 00 ) adds a unit margin above this root.

The following procedure computes A ̂ 10 ( a 00 ) (and thus the boundary of D) on a grid of a 00 values and can also be used to test membership of a given point ( a 00 , a 10 ) in D.

Algorithm 1 Computing A ̂ 10 ( a 00 ) and the boundary of D.

3.6.1 Reproducibility

Algorithm 1 is straightforward to implement with standard root bracketing. For interactive visualization of the boundary and the stability region, see the GeoGebra script mentioned in Section 5.

3.6.2 Numerical validation near the boundary

To validate Algorithm 1—and, indirectly, the characterization of Theorem 29—against an independent method, we computed the rightmost characteristic roots of Equation 4 using Chebyshev pseudospectral discretization of the infinitesimal generator (Breda et al., 2015) (N = 60 collocation points) at points placed just above and just below the computed envelope. Table 2 reports the results for ( a 01 , a 11 ) = ( 2,5 ) at two abscissae: a 00 = 5, where the fiber A 00 − 1 ( a 00 ) = { 2.7981 } is a singleton, and a 00 = 30, where it contains the three frequencies { 3.8969 , 6.6730 , 8.2233 } and the envelope is a genuine maximum over branches, attained at w * = 3.8969. In both cases, the rightmost root computed at the boundary point ( a 00 , A ̂ 10 ( a 00 ) ) is purely imaginary to machine precision, with imaginary part equal to the maximizing frequency w *, as predicted by Lemma 2; and its real part changes sign transversally across the boundary, decaying linearly with the distance ε, consistent with the nonvanishing crossing rate in Corollary 7. The numerical results in this study were obtained using Python 3 (NumPy) implementations of Algorithm 1 and of the pseudospectral discretization. All figures were generated using GeoGebra (International GeoGebra Institute, 2013) and edited using Inkscape, except the conceptual illustration of the supremum envelope (Figure 1), which was produced entirely using Inkscape.

TABLE 2. Validation of the computed boundary against pseudospectral root computations, for ( a 01 , a 11 ) = ( 2,5 ).

For each abscissa a 00, the table shows the real part of the rightmost characteristic root at the points ( a 00 , A ̂ 10 ( a 00 ) ± ε ). At the boundary point itself ( ε = 0 ), the computed rightmost root is Re s < 1 0 − 13 with Im s = w *.

4 Illustrative examples

Example 1 (Delayed position feedback without damping). Consider the particular case a 01 = a 10 = 0 (no instantaneous damping and no delayed position feedback). Then, A 00 ( w ) = − a 11 w ⁡ sin ( w ) + w 2 and A 10 ( w ) = − a 11 ⁡ cos ( w ). By Theorem 29 and Equation 12, the system is asymptotically stable if and only if 0 > A ̂ 10 ( a 00 ). This holds when 0 > A 10 ( w ) = − a 11 ⁡ cos ( w ), i.e., when w satisfies r k < w < r k + 1 / 2, where r k = π ( 4 k − sgn ( a 11 ) ) / 2 and r k + 1 / 2 = r k + π = π ( 4 k + 2 − sgn ( a 11 ) ) / 2, k ∈ Z +, and a 11 ≠ 0. Evaluating A 00 at the boundary points, we obtain the conditions Note that if | a 11 | ≥ π, then inequality (Equation 20) has no solution, yielding the necessary condition | a 11 | < π. Finally, since ( a 00 , 0 ) ∈ D ⊂ Ω, it must hold that a 00 + a 01 > 0. Together with a 01 = 0, this implies a 00 > 0. The same stability region is reported in 7, Proposition 2.5; this case was also analyzed by Ernest Bellman and Kenneth (1963). The resulting stability region is shown in Figure 3: it consists of a countable family of disjoint lobes stacked along the a 00-axis, one per admissible frequency band in Equation 20. The lobes alternate between positive and negative a 11 as k increases—reflecting the dependence of the bands on sgn ( a 11 )—and shrink linearly in width as | a 11 | increasing, vanishing at | a 11 | = π, which makes the necessary condition | a 11 | < π directly visible in the figure.

Example 2 (Air-to-fuel ratio dynamics in an internal combustion engine). Consider the mathematical model of the air-to-fuel ratio dynamics in a fully warm port-fuel-injected internal combustion engine presented by Orlov et al. (2009), where more detailed information about this system and model can also be found: Since only output measurement is available, and the delay τ may also be uncertain (e.g., during the identification process), it is natural to consider the output feedback control law v ( t ) = − k y ( t ). Substituting into Equations 21–23 yields the closed-loop equation presented in (Equation 24): Equation 24 is of the form Equation 2. After the time-scaling described in Section 2, the normalized system is obtained, as presented in Equation 25:The parameters in Equation 25 are given by Equation 26: The structure of Equation 26 clarifies the role of the gain k in the closed loop. In the model of Orlov et al. (2009), α and α 0 are the damping and stiffness coefficients of the intake/fuel-film dynamics, respectively, β weighs the derivative action of the delayed input, and τ is the cycle-and-transport delay between injection and measurement. The local channel ( a 00 , a 10 ) = ( τ 2 α 0 , τ α ) is fixed by the plant—normalized instantaneous stiffness and damping—and is not affected by the controller, whereas k scales the delayed channel ( a 01 , a 11 ) = k a 00 ( 1 , ℓ ) homothetically: increasing the gain rigidly amplifies both the delayed-position and the delayed-velocity feedback. The dimensionless ratio ℓ = β / ( τ α 0 ) measures the relative weight of the delayed-velocity action with respect to the delayed-position action; larger ℓ makes the delayed channel more derivative-like (and, as shown below, more aggressive per unit gain). The value a 00 = τ 2 α 0 = 0.102 used in Figure 4 corresponds to the engine operating point identified by Orlov et al. (2009). From Equation 7 and Equation 26, we obtain Equation 27: By Theorem 29 and the definition of D (Equation 12), asymptotic stability requires a 10 > A 10 ( w * ) for all w * such that A 00 ( w * ) = a 00. A sufficient condition is that a 10 > A 10 ( w ) for all w ∈ R. Sincestability is guaranteed if This condition is delay-independent: for any 0 < k < α / β, the closed-loop system is asymptotically stable regardless of the delay value τ > 0. This is particularly useful when the delay is uncertain or time-varying. Furthermore, if the system parameters are not precisely known but lower and upper bounds α min and β max are available, then 0 < k < α min / β max guarantees stability [note that k = 0 reduces to uncontrolled dynamics and the parametrization (Equation 26) is not valid in that case]. Two remarks quantify the conservativeness of the bound k < α / β. First, the condition is sufficient but not necessary: it is obtained from Theorem 29 through three successive relaxations—the fiber constraint A 00 ( w * ) = a 00 in Equation 11 is decreased, replacing a maximum over finitely many crossing frequencies by a supremum over all w > 0; the bracket in Equation 27 is bounded by sup w > 0 [ sin ⁡ w / w − ℓ ⁡ cos ⁡ w ] ≤ 1 + ℓ, which is never attained; and ( 1 + ℓ ) is finally replaced by ℓ. Second, the exact necessary-and-sufficient threshold is precisely the boundary curve A ( k ) introduced below: as Figure 4 shows, A ( k ) is even negative over an intermediate range of gains—where the delayed feedback contributes effective damping—so that every gain in that range is stabilizing for any positive plant damping a 10 = τ α > 0, whereas the delay-independent bound certifies only k < a 10 / ( 1 + ℓ ) a 00, which shrinks with the damping. The gap between the two criteria is largest for small a 00 = τ 2 α 0 (delay short with respect to the plant time scale), where the crossing frequencies selected by the fiber constraint are far from the maximizer of the bracket; conversely, the relaxation sup w [ sin ⁡ w / w − ℓ ⁡ cos ⁡ w ] ≤ 1 + ℓ becomes relatively tight for large ℓ since the supremum then approaches ℓ. When the delay is well identified, the exact boundary A ( k ) should, therefore, be preferred for tuning; the delay-independent bound remains the appropriate tool under delay uncertainty. If the parameters are known, one can obtain a tighter (delay-dependent) characterization of the stability region. From Equations 26, 27, definethat is, A ( k ) = A ̂ 10 ( a 00 ) as defined in Equation 11. Since both A 00 ( w ) and A 10 ( w ) do not depend on a 10 for a given value of k, the stability region in the ( k , a 10 )-plane is the set of points above the curve A ( k ) (see Figure 4).

Example 3 (Three-parameter stability visualization by slicing). Theorem 29 provides a complete description of the stability region in the ( a 00 , a 10 )-plane for fixed ( a 01 , a 11 ). In applications, however, one often wishes to understand how this region deforms as one of the delayed coefficients varies. A practical approach is to compute a family of planar sections and “stack” them to approximate a portion of the stability boundary in a higher-dimensional parameter space. For instance, fix a 11 and sample a 01 on a grid. For each value of a 01, Algorithm 1 yields the boundary curve a 10 = A ̂ 10 ( a 00 ) in the ( a 00 , a 10 )-plane. The union of these curves over the sampled a 01 values provides a triangulated surface that approximates the boundary of the stable set in ( a 00 , a 01 , a 10 ). Analogously, fixing a 01 and sweeping a 11 yields an approximation of the stability boundary in ( a 00 , a 11 , a 10 ). Planar intersections of such a surface recover stability regions in other two-parameter planes, which is useful for controller tuning and sensitivity studies (see Figure 5 for a conceptual illustration).

5 Discussion and conclusion

5.1 Summary

Theorem 29 provides a complete, constructive characterization of the asymptotic stability region in the parametric space for the general second-order single-delay linear system (Equation 3). The stability region D is defined explicitly via the supremum envelope A ̂ 10. Since D is an epigraph in the a 10-direction, it is connected and easily visualizable: its boundary is the graph a 10 = A ̂ 10 ( a 00 ), and D is the region above this curve inside Ω. In addition, Algorithm 1 in Section 3.6 provides a practical procedure to compute the boundary a 10 = A ̂ 10 ( a 00 ).

5.2 Relation to previous work

Such systems arise naturally when a second-order plantis controlled using state or output feedback subject to a delay. Table 1 summarizes previously reported stability results for particular sub-cases of Equation 3. Each of these results can be recovered from Theorem 29 by specializing the parameters accordingly.

x ̈ t + a 1 x ̇ t + a 0 x t = u t − 1

Note: For completeness, we mention the related neutral-type extension s 2 + a 10 s + a 00 + ( a 21 s 2 + a 11 s + a 01 ) e − s τ = 0 (with τ ≥ 0), whose stability regions are studied by Yang (1993). This case involves a delayed highest-derivative term and requires additional well-posedness/regularity assumptions; it is therefore outside the scope of the present work, which focuses on the retarded model (Equation 3).

Compared to the studies by Cahlon and Schmidt (2004) and Cahlon and Schmidt (2007b), where sufficient conditions are derived for specific sign configurations of the parameters, our result imposes no restrictions on the signs of a 00, a 10, a 01, or a 11 and provides conditions that are both necessary and sufficient.

The comparison with the recent PID-oriented literature is also instructive. The results of Ma and Chen (2019) and Chen et al. (2022) deliver, for a given low-order plant, the largest delay tolerated by PID control—exactly, in the second-order case considered by Chen et al. (2022)—while Appeltans et al. (2022) develop strongly stabilizing PID controllers for general delay systems; in all three cases, the object of interest is a margin along the delay axis for structured gains. Theorem 29 is complementary: for each fixed delayed channel ( a 01 , a 11 ), it returns the whole set of stabilizing local coefficients, from which delay margins are recovered as one-dimensional sections (Remark 11). Similarly, the MID methodology (Mazanti et al., 2021) assigns a prescribed dominant root, i.e., it selects a distinguished point (or curve) inside the stability region, whereas the present result describes the region itself, including its boundary—the natural object for the robustness quantification discussed below.

It is also instructive to compare with the stability lobe diagrams well known in the machining dynamics literature [see, e.g. (Tamás Insperger and Gábor Stépán, 2011; Altintaş and Budak, 1995)]. For the particular sub-case a 10 = a 11 = 0 (undamped, no delayed velocity feedback), the stability region described by Theorem 29 reduces to the classical stability lobe structure. Our result generalizes these diagrams to include arbitrary damping and delayed velocity terms.

5.3 Applications and robustness

The analytical characterization of D enables the definition of robustness indices by measuring the distance from a particular operating point ( a 00 , a 10 ) to the boundary of D. As demonstrated in Example 2, the characterization can also yield simple, delay-independent stability conditions for particular parameter structures, which are especially valuable when the delay is uncertain.

For practical use, the stability region can be visualized with the aid of interactive geometry software such as GeoGebra (International GeoGebra Institute, 2013), for which a script is available at https://ggbm.at/nrf8q485.

5.4 Limitations and future work

The scope and limitations of the present characterization deserve an explicit statement. First, the result concerns a single discrete delay: as made precise by the scaling map of Remark 11, the normalized delayed channel ( a 01 , a 11 ) = ( τ 2 a ̄ 01 , τ a ̄ 11 ) is coupled to τ, so the gain chart in the ( a 00 , a 10 )-plane is computed per delay value; delay sweeps require recomputing the envelope, although each evaluation is inexpensive (Algorithm 1). Second, the envelope A ̂ 10 need not be continuous (Remark 4): it may exhibit downward jumps at fold points of Γ, so robustness indices defined as distances to the boundary should be computed with respect to the closure of ∂ D. Third, on the algorithmic side, the sign-change bracketing in Algorithm 1 may miss solutions of A 00 ( w ) = a 00 of even multiplicity (tangential contact); in practice, this is remedied by refining the grid step h and monitoring local minima of | A 00 ( w ) − a 00 |, but a certified implementation would require interval arithmetic or the root-counting techniques of Vyhlídal and Zítek (2009).

The approach presented here also exploits the fact that, for second-order systems, the parametric curve ( A 00 ( w ) , A 10 ( w ) ) lies in a two-dimensional space matching the dimension of the parameter plane ( a 00 , a 10 ): the crossing locus is a curve in the same plane as the free parameters, so its upper envelope is well defined. For higher-order systems, the stability boundary becomes a manifold of higher co-dimension, and the supremum envelope construction does not directly generalize without further development. Extending the present approach to third-order or higher odd-order systems (Cahlon and Schmidt, 2007b) is an interesting direction for future research, potentially in combination with computational techniques from Vyhlídal and Zítek (2009), Michiels and Niculescu (2014), Breda et al. (2015).

Another natural extension concerns systems with multiple or distributed delays, where the parametric structure of the stability boundary is considerably richer; (see Gu et al., 2005; Sipahi et al., 2011; Qi et al., 2019 for relevant frameworks). Finally, many engineered systems are additionally subject to impulsive effects or abrupt state jumps. The stability theory of impulsive time-delay systems is developed in a hybrid-systems framework [see, e.g. (Liu and Zhang, 2020)], in which the characteristic quasi-polynomial—the very object underlying the D-partition construction used here—is no longer available; whether an envelope-type description of stabilizing coefficient sets exists for impulsive delay systems is an open question beyond the scope of the present frequency-domain approach.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. 1. Mathematical review and proof refinement. All lemmas, propositions, corollaries, and their proofs were submitted iteratively to Claude acting as an adversarial mathematical reviewer. 2. Rewriting definitions, lemma statements, and proofs for clarity and elegance, applying consistent notation throughout. All the results, and the insights for their proofs, were provided by the authors.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Vinícius Piccirillo, Federal University of Technology, Brazil

Reviewed by: Denys Baranovskyi, Rzeszów University of Technology, Poland

Earige Tanveer, Kohsar University Murree, Pakistan