# On the Sensitivity of Characteristic Roots of a Class of Parameterized Delay-Differential Neutral Systems

César Fernando Méndez Barrios, Silviu-Iulian Niculescu, Alejandro

Martínez-González

► To cite this version:

César Fernando Méndez Barrios, Silviu-Iulian Niculescu, Alejandro Martínez-González. On the Sensitivity of Characteristic Roots of a Class of Parameterized Delay-Differential Neutral Systems. CoDIT 2024 - 10th International Conference on Control, Decision and Information Technologies, Jul 2024, La Valette, Malta. 10.1109/codit62066.2024.10708501. hal-04633786

HAL Id: hal-04633786

https://hal.science/hal-04633786v1

Submitted on 3 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the Sensitivity of Characteristic Roots of a Class of Parameterized Delay-Differential Neutral Systems

César-Fernando Méndez-Barrios<fnref n="1" />,*, Silviu-Iulian Niculescu<fnref n="2" />, and Alejandro Martínez-González<fnref n="3" />

**Abstract**— This paper focuses on the characterization of the asymptotic behavior of the critical characteristic roots for a class of delay-differential dynamical systems of neutral type whose coefficients smoothly depend on certain parameters. Such systems can be described by coupled delay-differential and delay-difference equations, and model time heterogeneity in propagation and transport phenomena. The asymptotic behavior of the characteristic roots is addressed by expressing the solutions as a convergent Puiseux series, which facilitates handling multiple solutions. Particular attention is paid to the way the parameters affect the stability of the delay-difference operator. Illustrative examples complete the presentation and show the effectiveness of the proposed method.

## I. INTRODUCTION

Propagation and transport are typical phenomena whose mathematical models can be described by delay-differential equations (DDEs) of neutral type (see, e.g., [1], [2], [3] and [4], [5] and the references therein). For further examples, we refer to [6]–[7]. As pointed out in [2] and [3], the exponential stability of the null solution of the associated delay-difference equation (in continuous time) represents a necessary condition for the exponential stability of the (null solution of the) corresponding DDEs of neutral type (see also [1]). In the linear case, except for the point spectrum, we will also have an essential spectrum that cannot be removed with a bounded perturbation. Thus, in the case of neutral DDEs, the stability analysis is more involved, and a deeper comprehension is required to understand the sensitivity of the spectrum with respect to the system parameters.

Motivated by the previous observations, this paper addresses the problem above for parameterized linear DDEs of neutral type, and, more precisely, we are interested in the behavior of critical characteristic roots of finite multiplicity and related splitting characterization. Under appropriate assumptions, this paper extends some of the ideas presented in [8] in the retarded case to deal with multiple roots located on the imaginary axis as a function of the system parameters (including the delay). For the computation of the characteristic roots located on the imaginary axis, we use the ideas based on matrix pencil properties proposed by [9] (see also [10] for further discussions). For a different

approach and related methodology, we refer to [11]–[12] and the references therein.

Our approach makes use of the Puiseux Theorem and, inspired by ideas developed in [13], we extend such ideas to a more general class of dynamical systems. More precisely, we compute the leading terms of the solution roots (represented in the form of Puiseux series) around the critical multiple roots. As expected, particular attention will be paid to the spectral properties of the associated delay-difference equation as a function of the parameters that appears in the neutral case but not in the retarded case. To the best of the authors' knowledge, there do not exist similar results in the open literature, and it represents a novelty.

Throughout the paper, the following standard notations will be adopted: $\mathbb{C}$ is the set of complex numbers, $i := \sqrt{-1}$, for $z \in \mathbb{C}$, $\text{Arg}\{z\} \in (-\pi, \pi]$ denotes the main argument. The unitary open (closed) disk will be denoted by $\mathbb{D}$ ($\mathbb{D}$). For a matrix $A$, its spectrum is denoted by $\sigma(A)$, and the $i$-th eigenvalue by $\lambda_i(A)$. If $P$ is a polynomial (quasi-polynomial), then $\sigma(P)$ denotes the set of roots of $P$, and its degree is denoted by $\text{deg}(P) = n, n \in \mathbb{N}$.

## II. PREREQUISITES AND PROBLEM FORMULATION

### A. Neutral Time-Delay Systems

Under appropriate initial conditions, consider the Linear Time-Invariant Delay Systems with a single discrete delay $\tau > 0$ of neutral type:

$$ \dot{y}(t) - B_d(p)\dot{y}(t - \tau) = A(p)y(t) + A_d(p)y(t - \tau). \quad (1) $$

The above representation is inspired by [14], in (1) the matrix functions $A, A_d, B_d: \mathcal{P} \mapsto \mathbb{R}^{n \times n}$ are assumed to be continuously differentiable functions of $p$, where $\mathcal{P} \subset \mathbb{R}$ is an appropriate open set. The characteristic function of (1) is denoted by $f: \mathbb{C} \times \mathbb{R}_+ \times \mathcal{P} \to \mathbb{C}$ and is defined as

$$ f(s; \tau, p) := \sum_{k=0}^{q} p_k(s; p) e^{-k\tau s} \quad (2) $$

where $p_k(s; p)$ are polynomial functions of $s$ given by

$$ p_0(s; p) \triangleq s^n + \sum_{i=0}^{n-1} a_{0i}(p) s^i, \quad p_k(s; p) \triangleq b_k(p) s^n + \sum_{\ell=0}^{n-1} a_{k\ell}(p) s^\ell, $$

and the coefficients $a_{k\ell}$ and $b_k$ are assumed to be continuously differentiable functions of $p \in \mathcal{P}$. It is well known that in order to achieve the asymptotic stability of (1) for arbitrary delay values, it is necessary to preserve the stability of the

<footnote>*Corresponding author.</footnote>

<footnote n="1">¹C.-F. Méndez-Barrios is with Universidad Autónoma de San Luis Potosí (UASLP), Facultad de Ingeniería, Dr. Manuel Nava No.8, cerfranfer@gmail.com</footnote>

<footnote n="2">²S.-I. Niculescu is with Université Paris-Saclay, CNRS, CentraleSupélec, Inria (DISCO team), Laboratoire des Signaux et Systèmes (L2S, UMR CNRS 8506) F-91192, Gif-sur-Yvette, France.</footnote>

<footnote n="3">³A. Martínez-González is with IPICyT, Mexico.</footnote>
neutral chain (see, for instance [9]). Such a neutral part is characterized by the following difference equation:

$$y(t) + \sum_{k=1}^{q} b_k(p)y(t - k\tau) = 0. \quad (3)$$

In this vein, the stability analysis of (3) can be performed by means of the following result:

**Proposition 2.1 ([9]):** Let $\mathcal{P}$ be an open and non-empty subset of $\mathbb{R}$. Then, the difference equation (3) will be asymptotically stable for all $\tau \ge 0$ and all $p \in \mathcal{P}$ if and only if the function $N$ defined by

$$N(s; p) = 1 + \sum_{k=1}^{q} b_k(p) e^{-k\tau s}, \quad (4)$$

has all its roots in $\mathbb{C}_-$.

Bearing in mind the previous result, let $p^* \in \mathcal{P}$ and $\tau^* > 0$ be fixed values. Then, (1) will be asymptotically stable if Proposition 2.1 holds and all the solutions of quasi-polynomial (2) are located at $\mathbb{C}_-$ (see, for instance, [9]).

## B. Asymptotic Behavior

As mentioned, the behavior of the solutions on the imaginary axis for smooth variations of the parameter $p$ (or $\tau$) becomes of core importance. Furthermore, the task of describing such behavior gets much more complex when the solution is multiple ($m \in \mathbb{N}$, s.t. $m > 1$). A key notion to describe such behavior is to express the solution as Puiseux series [15]. In this regard, by denoting the parameter $p$ or $\tau$ by $\xi$, the equation $f(s; \xi) = 0$ defines a solution curve $C \in \mathbb{C}^2$ which is composed by the finite union of $r$-branches $s_j(\xi^{1/m_j})$, each of these branches can be expressed as a Puiseux series:

$$s_{j,\sigma}(\xi) = c_{j,\sigma}\xi^{\frac{1}{m_j}} + o\left(|\xi|^{\frac{1}{m_j}}\right), \quad j = 1, \dots, r, \quad \sigma = 1, \dots, m_j,$$

where each branch has multiplicity $m_j$, such that $m = m_1 + m_2 + \dots + m_r$. In the case when $r = 1$, then $s_{j,\sigma}$ and $c_{j,\sigma}$ will be simply denoted by $s_\sigma$ and $c_\sigma$, respectively.

## C. Problem Formulation

As mentioned in the Introduction, this note focuses explicitly on the following problems:

(i) first, for a given quasi-polynomial $f(s; \tau, p)$ and a known simple solution $(i\omega^*, \tau^*, p^*) \in \mathbb{C} \times \mathbb{R}_+ \times \mathbb{R}$ find all coefficients $\gamma_j \in \mathbb{C}$, $j \in \mathbb{N}$ of the power series expansion of the solution $s(p)$, i.e,

$$s(\xi) = i\omega^* + \sum_{j=1}^{\infty} \gamma_j(\xi - \xi^*)^j,$$

(ii) second, for a known multiple solution $(i\omega^*, \tau^*, p^*) \in \mathbb{C} \times \mathbb{R}_+ \times \mathbb{R}$, find the first coefficients of the Puiseux series expansion of the solution $s(\xi)$,

$$s(\xi) = i\omega^* + \gamma_1(\xi - \xi^*)^{\frac{1}{m}} + o\left(|\xi - \xi^*|^{\frac{1}{m}}\right),$$

and,

(iii) finally, find the stability crossing directions, that is, determine whether the solution $s(\xi)$ enters to the right half-plane of the complex plane (RHP), or to the left half-plane of the complex plane (LHP) for $\xi > \xi^*$,

where $\xi$ denotes one of the parameters $\tau$ or $p$, according to the context. Observe that in the case of Problem-(i), our aim is to give a complete characterization of the solution $s(\xi)$.

# III. MAIN RESULTS

It is well known that, even though the characteristic roots may all have negative real parts, it is still possible for some solutions to be unbounded (see, for instance, [9]). Keeping in mind such a fact, we will consider conditions ensuring that the characteristic function has continuous bounded roots. In this manner, we will take a systematic stability analysis concerning the behavior of the Critical Roots located on the imaginary axis of the complex plane, either simple or multiple, when the parameter $p$ (or the delay $\tau$) is under small smooth variations.

## A. Critical Values Characterization

It is well known that when (3) is stable, the continuity of the solutions $f(s; \tau, p) = 0$ with respect to $p$ or $\tau$ is ensured. This property has an important implication, that is, the change in stability with respect to the delay value is determined by the asymptotic behavior of the critical roots $s = i\omega$. Hence, it is of vital importance to determine conditions on $\tau$ or $p$ that allow detecting critical solutions, that is, solutions on the imaginary axis, or in other words, to find $\tau$ of $p$ for which there exist $s^* = i\omega$ such that $f(i\omega; \tau, p) = 0$. The following result allows this procedure.

**Proposition 3.1 ([9]):** For a fixed $p^* \in \mathbb{R}$ assume that conditions in Proposition 2.1 holds. For $k = 0, 1, \dots, n-1$, introduce

$$T_n(p^*) := \begin{bmatrix} 1 & 0 & \cdots & 0 \\ b_1(p^*) & 1 & \cdots & 0 \\ \vdots & \ddots & \ddots & \vdots \\ b_{q-1}(p^*) & b_{q-2}(p^*) & \cdots & 1 \end{bmatrix},$$

$$T_k(p^*) := \begin{bmatrix} a_{0i}(p^*) & 0 & \cdots & 0 \\ a_{1i}(p^*) & a_{0i}(p^*) & \cdots & 0 \\ \vdots & \ddots & \ddots & \vdots \\ a_{q-1,i}(p^*) & a_{q-2,i}(p^*) & \cdots & a_{0i}(p^*) \end{bmatrix},$$
$$
H_n(p^*) := \begin{bmatrix} b_q(p^*) & b_{q-1}(p^*) & \cdots & b_1(p^*) \\ 0 & b_q(p^*) & \cdots & b_2(p^*) \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & b_q(p^*) \end{bmatrix},
$$

$$
H_k(p^*) := \begin{bmatrix} a_{qi}(p^*) & a_{q-1i}(p^*) & \cdots & a_{1i}(p^*) \\ 0 & a_{qi}(p^*) & \cdots & a_{2i}(p^*) \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_q(p^*) \end{bmatrix},
$$

$$
P_k(p^*) := \begin{bmatrix} (\mathbf{i})^i T_i(p^*) & (\mathbf{i})^i H_i(p^*) \\ (-\mathbf{i})^i H_i^T(p^*) & (-\mathbf{i})^i T_i^T(p^*) \end{bmatrix}, k = 0, \dots, n.
$$

Define further,

$$
P(p^*) := \begin{bmatrix} 0 & I & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & I \\ -P_n^{-1}P_0 & -P_n^{-1}P_1 & \cdots & -P_n^{-1}P_{n-1} \end{bmatrix}.
$$

$$
F(s) := \begin{bmatrix} 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \\ -p_0(s; p) & -p_1(s; p) & \cdots & -p_{q-1}(s; p) \end{bmatrix} \quad (5)
$$

$$
G(s) := \mathrm{diag}(1, \dots, 1, p_q(s; p)). \quad (6)
$$

Then, the critical delay values are given by

$$
\tau_k := \frac{\alpha_{k,i}}{\omega_k}, \quad i = 1, \dots, m, \quad k = 1, \dots, l,
$$

where $m \le q$, $l \le 2nq$, $0 \ne \omega_k \in \sigma(P) \cap \mathbb{R}_+$ and $\alpha_{k,i} \in [0, 2\pi)$ such that satisfy $e^{-j\alpha_k} \in \sigma(F(j\omega_k), G(j\omega_k))$.

**Remark 3.1:** It is worth mentioning that Proposition 3.1 results as a slight modification to *Theorem 1* presented in [16].

According to Proposition 3.1, we need to provide the stability
of the neutral part, that is, guaranteeing the stability of N.
In this regard, our aim is to determine the set of all critical
values

$$
\mathcal{P} := \{ p_j^* \in \mathbb{R} : j \in \mathbb{N} \text{ and } N(i\omega, p_j^*) = 0 \text{ for some } \omega \ge 0 \}.
$$

Since $\mathcal{P}$ is a countable set, without any loss of generality,
for $\mathcal{P} \neq \emptyset$, we will assume that $\mathcal{P}$ is an ordered set, that
is, if $p_j^*, p_k^* \in \mathcal{P}$ with $j < k$, then $p_j^* < p_k^*$. Under these
considerations, if $\mathcal{P} \neq \emptyset$ and, if for some $p^* \in \mathbb{R}$ we have
that $N$ is stable, then, $N$ will remain stable for any $p \in$
$(p_j, p_{j+1})$, where $p_j$ satisfy $p_j < p^* < p_{j+1}$. To determine
the set $\mathcal{P}$, let us introduce the following parameter-dependent
matrix $W_N : \mathbb{R} \to \mathbb{R}^{q \times q}$:

$$
W_N(p) := T_n(p) T_n^T(p) - H_n^T(p) H_n(p) . \quad (7)
$$

We have the following :

*Proposition 3.2:* Consider the neutral quasi-polynomial
(4). The following statement holds:

$p^*$ is a critical point if and only if

$$
(i) \ 0 \in \sigma(W_{\hat{N}}(p^*));
$$

$$
(ii) \sigma(\hat{N}(z; p^*)) \cap \partial\mathbb{D} \neq \emptyset;
$$

where the function $\hat{N}(z; p)$ is a polynomial defined by

$$
\hat{N}(z; p) := z^q + \sum_{k=1}^{q} b_k(p) z^{q-k}. \quad (8)
$$

The stability of the characteristic function N(s, p) plays a fundamental role in our analysis because it allows the continuity of the roots of f(s; τ, p). In this situation, with N(s, p) stable, the crossing directions of critical roots allows the stability analysis of the neutral quasi-polynomial f.

B. Asymptotic Characterization

Let us consider the asymptotic behavior of simple and
multiple critical roots for both f(s; τ, p) and N̂(z; p). With
the aim of presenting a unified approach, we will denote by η
to refer to either f or N̂; similarly, μ will indicate a solution
on D or ∂C+, and, according to the context, ξ will represent
either p or τ.

1) *Simple Solutions*: as a first step, let us derive the complete characterization of a simple solution μ(ξ). To this end, the following result facilitates this task.

*Proposition 3.3:* Let $\mu^* = \mu_0 \in \mathbb{C}$ at $\xi = \xi^*$ be a simple critical solution of the function $\eta(\mu, \xi)$. Then, the complete description of $\mu(\xi)$ around the critical point is given by

$$
\mu(\xi) = \mu_0 + \sum_{k=1}^{\infty} c_k (\xi - \xi^*)^k, \quad (9)
$$

where the coefficient $c_k$ are being computed as follows:

$$
c_k = - \left( \frac{\partial \eta}{\partial \mu} \right)^{-1} \left[ \frac{\partial^k \eta}{\partial \xi^k} + \sum_{\ell=1}^{k-1} \frac{\partial^k \eta}{\partial \xi^{k-\ell} \partial \mu^\ell} c_\ell + \right. \\
\left. \qquad + \sum \frac{j!}{j_1! j_2! \cdots j_r!} \frac{\partial^{\ell+j} \eta}{\partial \mu^\ell \partial \xi^j} c_1^{j_1} c_2^{j_2} \cdots c_r^{j_r} \right], \quad (10)
$$

where the summation take values over

$$
\begin{align*}
& j \in \{2, 3, \dots\}, \quad r \in \{1, 2, \dots\}, \\
& j_1, j_2, \dots, j_r \in \{0, 1, \dots\}, \quad \ell \in \{0, 1, \dots\}
\end{align*}
$$

such that $\ell + j_1 + 2j_2 + \cdots + rj_r = k$.

*Remark 3.2:* It can be seen from Proposition 3.3 that the higher-order terms are determined by following a recursive procedure.

*Remark 3.3:* As we will see later, the complete descrip-
tion of the asymptotic behavior of the solutions will allow
a better comprehension of the degenerate cases. Specifically,
such a description will be very useful in determining the
crossing directions.

2) *Multiple Solutions*: let us consider now the asymptotic behavior of a *m*-multiple critical root $\mu^*$ of $\eta$ at $\xi^*$. In this case, it is possible to characterize the root locus of $\eta$ around a multiple root by its branches (see, for instance, [15]), where each of these branches can be expressed as a Puiseux series:

$$
\mu_{j,\sigma}(\xi) = c_{j,\sigma} \xi^{\frac{1}{m_j}} + o\left(|\xi|^{\frac{1}{m_j}}\right), \quad (11)
$$
with $j = 1, \dots, r$, $\sigma = 1, \dots, m_j$. The parametrization of each branch has (partial) multiplicity $m_j$, such that $m = m_1 + m_2 + \dots + m_r$. In the case when $r = 1$, then $\mu_{j\sigma}$ and $c_{j,\sigma}$ will be simply denoted by $\mu$ and $c_\sigma$, respectively.

Now, to derive a characterization of a multiple root around $\mu^*$, we consider the structure (11). For such an end, let us introduce the constant $\rho \in \mathbb{N}$, which has a similar to the multiplicity $m$, that is, let $\rho$ be the natural number satisfying:

$$ \left. \frac{\partial \eta}{\partial \xi} \right|_{(\mu^*; \xi^*)} = \dots = \left. \frac{\partial^{\rho-1} \eta}{\partial \xi^{\rho-1}} \right|_{(\mu^*; \xi^*)} = 0, \quad \left. \frac{\partial^\rho \eta}{\partial \xi^\rho} \right|_{(\mu^*; \xi^*)} \neq 0. \quad (12) $$

Bearing the above notation in mind, the following result allows deriving the first terms of the Puiseux series.

**Proposition 3.4:** Let $\mu = \mu^*$ be a $m$-multiple root of $\eta(\mu; \xi)$ at $\xi = \xi^*$. Then, the $m$-zeros can be expanded as

$$ \mu_{j,\sigma}(\xi) = \mu^* + c_{j,\sigma} (\xi - \xi^*)^{\beta_j} + o(|\xi - \xi^*|^{\beta_j}), $$

for $j = 1, 2, \dots, r$, $\sigma = 1, \dots, m_j$ and $m = m_1 + \dots + m_r$. The characterization of the splitting behavior is given by the following cases:

(i) if $\rho = 1$, then there exists only one ramification, that is, $r = 1$, its exponent is $\beta_1 = 1/m$ and, its corresponding coefficient $c_\sigma$ is given by:

$$ c_{\sigma} = \left| m! \left( \frac{\frac{\partial \eta}{\partial \xi}}{\frac{\partial^m \eta}{\partial \mu^m}} \right)^{\frac{1}{m}} e^{i \frac{(2\sigma-1)\pi + \theta}{m}} \right|, \quad (13) $$

$$ \text{where } \theta := \operatorname{Arg} \left\{ \frac{\partial \eta}{\partial \xi} / \frac{\partial^m \eta}{\partial \mu^m} \right\} \text{ and } \sigma = 1, 2, \dots, m. $$

(ii) if $\frac{\partial^2 \eta}{\partial \mu \partial \xi} \neq 0, 1 < \rho < \infty$ and $m \ge 2$, the first ramification posses exponent $\beta_1 = 1/m_1$ and coefficient:

$$ c_{1,\sigma} = \left| m_1! \frac{\frac{\partial^2 \eta}{\partial \xi \partial \mu}}{\frac{\partial^m 1 \eta}{\partial \mu^m}} \right|^{\frac{1}{m_1}} e^{i \frac{(2\sigma-1)\pi + \theta}{m_1}}, \quad (14) $$

where $m_1 := m - 1$, $\theta := \operatorname{Arg} \left\{ \frac{\partial^2 \eta}{\partial \xi \partial \mu} / \frac{\partial^m 1 \eta}{\partial \mu^m} \right\}$ and $\sigma = 1, \dots, m_1$. For the second ramification, the corresponding coefficients are given by:

$$ c_{2,1} = - \frac{\frac{\partial^\rho \eta}{\partial \xi^\rho}}{\frac{\partial^2 \eta}{\partial \xi \partial \mu}}, \quad \text{and} \quad \beta_2 = \rho - 1, \quad m_2 = 1. \quad (15) $$

Otherwise $\mu_{2,1}$ is an invariant root, i.e., $\mu_{2,1} = \mu^*$ for all $\xi$.

(iii) if $\rho = m$ and $\frac{\partial \eta}{\partial \xi^\alpha \partial \mu^\beta} = 0$ for some $\alpha, \beta \in \mathbb{N}$ such that $\alpha + \beta > m$, then there is only one ramification with exponent $\beta_1 = 1$, where its coefficients are given by:

$$ c_{\sigma} = \left| m! \left( \frac{\frac{\partial^m \eta}{\partial \xi^m}}{\frac{\partial^m \eta}{\partial \mu^m}} \right)^{\frac{1}{m}} e^{i \frac{(2\sigma-1)\pi + \theta}{m}} \right|, \quad (16) $$

$$ \text{where } \theta := \operatorname{Arg} \left\{ \frac{\partial^m \eta}{\partial \xi^m} / \frac{\partial^m \eta}{\partial \mu^m} \right\} \text{ and } \sigma = 1, 2, \dots, m. $$

### C. Crossing Directions Characterization

To derive a complete characterization of the stability properties of the critical solutions, it will be necessary to understand its local behavior. In this context, the remaining part of this section will analyze the crossing directions of a critical solution as the parameter $\xi$ crosses the critical value $\xi^*$ in the increasing direction. In the seek for simplicity, the analysis of the neutral part will be carried out by means of the polynomial $\tilde{N}$ introduced in (8), which implies that we will be looking for the behavior of the solution $z(p)$ of the equation

$$ \tilde{N}(z; p) = 0, \quad (17) $$

on the unit circle. In this regard, we will examine such a solution to determine when it crosses the unit circle in the inner or outer direction or even when $z(p)$ touches it tangentially. In summary, the conduct of $z(p)$ around $z(p^*)$ will have one of the following behaviors:

(i) it crosses the unit circle;

(ii) it is tangent to a point on the unit circle;

(iii) it intersects the unit circle with other solutions;

Some of the above possibilities are described in Fig. 1. It is worth to mention that similar behaviors may occur for the solutions on the imaginary axis.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Some possible behaviors of $z(p)$. (Upper) cases (i)-(ii). (Lower) case (iii).</figcaption></figure>

1) Simple Solutions: the following result seeks to describe such behavior in such a way that crossing direction is described.

**Proposition 3.5:** Under the assumption that the critical solutions of $\tilde{N}$ are simple, the following statements are equivalent:

(i) The root $z^* = e^{i\theta^*}$ is crossing $\partial\mathbb{D}$ towards instability (stability).
(ii) The following inequality holds:

$$ \left. \frac{d}{dp} \left| z \right| \right|_{p=p^*} > 0 \quad (< 0), $$

for any $p$ sufficiently close to $p^*$, but $p > p^*$.

(iii) The following inequality holds:

$$ \Re \left\{ \frac{d \hat{N}(e^{i\theta^*}; p)}{dp} \bigg|_{p=p^*} - \frac{d \hat{N}(z; p^*)}{dz} \bigg|_{z=e^{i\theta^*}} \right\} < 0 \quad (> 0), \quad (18) $$

then for any $p$ sufficiently close to $p^*$ but $p > p^*$.

It is worth mentioning that under some situations (18) vanishes, in such a case, the following result allows for providing a second-order analysis.

**Proposition 3.6:** Let $p = p^*$ be a critical parameter, such that $z^* = e^{i\theta^*}$ corresponds to the critical solution of $\hat{N}$. Under the assumption that $z^*$ is simple, the following statements are equivalent:

(i) The solution $z(p)$ stays outside (inside) of the unit circle $\partial\mathbb{D}$.

(ii) The following inequality holds:

$$ \left. \frac{d^2 |z|}{dp^2} \right|_{p=p^*} > 0, \quad (< 0) \quad (19) $$

for any $p$ sufficiently close to $p^*$, but $p > p^*$.

(iii) The following inequality holds:

$$ \Re \left\{ \frac{\varphi(z; p)}{z \left( \frac{\partial \hat{N}(z; p^*)}{\partial z} \right)^3} \bigg|_{\substack{z=z^* \\ p=p^*}} \right\} + \Im \left\{ \frac{\frac{\partial \hat{N}(z; p)}{\partial p} \bigg|_{p=p^*}}{z \frac{\partial \hat{N}(z; p^*)}{\partial z} \bigg|_{z=z^*}} \right\} > 0, (< 0) $$

where

$$ \varphi(z; p) := 2 \frac{\partial^2 \hat{N}}{\partial z \partial p} \frac{\partial \hat{N}}{\partial p} \frac{\partial \hat{N}}{\partial z} - \frac{\partial^2 \hat{N}}{\partial z^2} \left( \frac{\partial \hat{N}}{\partial p} \right)^2, $$

for any $p$ sufficiently close to $p^*$ but $p > p^*$.

To give a complete description of the asymptotic behavior of simple roots, the following result gives a general framework.

**Proposition 3.7:** Let $\xi = \xi^*$ be a critical parameter and $\mu^* \in \mathbb{C}$ its corresponding critical solution, such that $\Re(c_j) \equiv 0$ or $\cos(\Arg\{c_j\} - \Arg\{\mu^*\}) \equiv 0$, for $j = 1, 2, \dots, k-1$, where $c_j$ are the coefficients given in (9). Then, the root $\mu$ for sufficiently close to $\xi$ is crossing towards instability (stability) if

$$ (i) \text{ if } \eta = f: \qquad \Re(c_k) > 0 \quad (< 0). \quad (20) $$

$$ (ii) \text{ if } \eta = \hat{N}: \qquad \cos(\text{Arg}\{c_k\} - \text{Arg}\{\mu^*\}) > 0 \quad (< 0). \quad (21) $$

for any $\xi$ sufficiently close to $\xi^*$ but $\xi > \xi^*$.

2) **Multiple Solutions:** finally, and to complete the analysis, the following results give the characterization of the crossing direction of MCR's.

**Proposition 3.8:** Let $\xi^*$ be a $m$-multiple root of $\eta(\mu, \xi)$. For $\xi > \xi^*$ sufficiently close $\xi^*$, the characteristic root $\mu_{j,\sigma}$ (3.4) will enter the unstable region (stable region) if

$$ \cos(\text{Arg}\{c_{j\sigma}\} - \theta^*) > 0 (< 0), \quad \sigma = 0, 1, \dots, m_j, \quad (22) $$
where $\theta^* = \text{Arg}\{\mu^*\}$.

## IV. ILLUSTRATIVE EXAMPLES

In order to illustrate the effectiveness of the proposed approach, several numerical examples are proposed.

**Example 4.1:** As a first example, let us consider the quasi-polynomial $f(s; \tau, p) :=$

$$ p_0(s) + p_1(s, p)e^{-\tau s} + p_2(s, p)e^{-2\tau s} + p_3(s, p)e^{-3\tau s} \quad (23) $$

where

$$ \begin{align*} p_0(s) &:= s^4 + 2s^3 + 5s^2 + 3s + 2, \\ p_1(s, p) &:= (p + \frac{3}{5}) s^4 + s^2 + 2, \\ p_2(s, p) &:= \left(\frac{99p}{50} + \frac{11}{100}\right) s^4 + s^2 + s + 2, \\ p_3(s, p) &:= \left(\frac{47p^2}{125} + \frac{3}{500}\right) s^4 + 2s^3 + 5s. \end{align*} $$

From (23) we have that $\hat{N}$ is given by:

$$ \hat{N}(z; p) := z^3 + (p + \frac{3}{5})z^2 + \left(\frac{99p}{50} + \frac{11}{100}\right)z + \left(\frac{47p^2}{125} + \frac{3}{500}\right). \quad (24) $$

Thus, as a first step, we determine the stability of the neutral part. To such an end, we apply Proposition 3.2 and Proposition 3.5. Table I summarizes the results.

<table><caption>TABLE I: Critical parameters for the delay-difference operator in quasi-polynomial (23).</caption><thead><tr><th>p<sub>k</sub></th><th>Proposition 3.2-(ii)</th><th>{z ∈ C | N̂(z; p*) = 0}</th><th>Re {<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">dN̂/dp</div><div style="padding-top: 1px;">z | dN̂/dz</div></div> |<sub>z=z*</sub></div></th></tr></thead><tbody><tr><td>p<sub>1</sub></td><td>-7.300382509</td><td>{-2.45723, 1, 8.15761}</td><td>0.101427 (+)</td></tr><tr><td>p<sub>2</sub></td><td>-0.6251494057</td><td>{-1.11235, 0.137498, 1}</td><td>1.377616 (+)</td></tr><tr><td>p<sub>3</sub></td><td>-0.4400048822</td><td>{-1, 0.107582, 0.732423}</td><td>0.683178 (+)</td></tr><tr><td>p<sub>4</sub></td><td>0.5</td><td>{-0.5 ± i0.866025, -0.1}</td><td>-0.847032 (-)</td></tr><tr><td>p<sub>5</sub></td><td>3.046387861</td><td>{-1.32319 ± i1.32084, -1}</td><td>-0.708941 (-)</td></tr></tbody></table>

From the above, it follows that the set $\mathcal{P}$ is given by

$$ \mathcal{P} := \{-7.300382, -0.625149, -0.440004, 0.5, 3.046387\}. $$

Moreover, since $\sigma\{\hat{N}(z; 10)\} = \{-2.747, 1.256, 10.891\}$, Table I allows concluding that $I_s = (-0.4400048, 0.5)$ corresponds to the stable interval for $\hat{N}$. Figure 2 illustrates the behavior of $\sigma(\hat{N}(z; p))$ for $p \in (-\frac{3}{5}, \frac{3}{4})$. Now, by taking $p^* = 0 \in I_s$ we apply Propositions 3.1, 3.3 and 3.7 to obtain the results summarized in the Table II. Since $f(0, 0, 0)$ is stable, Table II allows us to conclude that system (23) will be asymptotically stable for any $\tau \in [0, 0.135612)$.

**Example 4.2:** As a final example, consider the following quasi-polynomial:

$$ f(s; \tau, p) := (s + a(p)) + \left(\frac{1}{5}s + b(p)\right)e^{-\tau s} + \left(\frac{2}{5}s + d(p)\right)e^{-3s\tau} $$
<figure><img src="image_3.png" /><figcaption>Fig. 2: Behavior of $z(p)$ for $p \in (-\frac{3}{5}, \frac{3}{4})$.</figcaption></figure>

<table><caption>TABLE II: Critical delays of quasi-polynomial (23).</caption><thead><tr><th colspan="2">Proposition 3.1</th><th colspan="2">Proposition 3.3</th></tr><tr><th><i>i</i>&omega;<sub>k</sub></th><th>&tau;<sub>k</sub></th><th>c<sub>1</sub></th><th>Proposition 3.7</th></tr></thead><tbody><tr><td>&plusmn;<i>i</i>3.017122</td><td>1.580251267</td><td>0.135678 &plusmn; <i>i</i>1.690265</td><td>0.1356 (+)</td></tr><tr><td>&plusmn;<i>i</i>1.872312</td><td>2.868857313</td><td>-0.081335 &plusmn; <i>i</i>0.441718</td><td>-0.0813 (-)</td></tr><tr><td>&plusmn;<i>i</i>1.443155</td><td>0.135612412</td><td>0.519451 &plusmn; <i>i</i>0.489140</td><td>0.5194 (+)</td></tr><tr><td>&plusmn;<i>i</i>0.869184</td><td>3.340121652</td><td>0.032804 &plusmn; <i>i</i>0.206203</td><td>0.0328 (+)</td></tr><tr><td>&plusmn;<i>i</i>0.273967</td><td>13.12182504</td><td>-0.000922 &plusmn; <i>i</i>0.018632</td><td>-0.0009 (-)</td></tr></tbody></table>

where the coefficients $a, b, d$ are given by:

$$a(p) = -\frac{1}{8} \csc\left(\frac{3p}{8}\right)\left(8b \sin\left(\frac{p}{4}\right) + 5\left(\cos\left(\frac{p}{4}\right) + 5\cos\left(\frac{3p}{8}\right) + 2\right)\right),$$

$$b(p) = \frac{1}{8} \left( 25 \csc\left(\frac{p}{4}\right) \left( 8 \sin\left(\frac{3p}{8}\right) - \frac{p}{5} \left( \cos\left(\frac{p}{4}\right) + 6 \right) \right) + 40 \right),$$

$$d(p) = \frac{1}{8} \csc\left(\frac{3p}{8}\right) \left( 5\left(\cos\left(\frac{p}{8}\right) + 2\cos\left(\frac{3p}{8}\right) + 5\right) - 8b \sin\left(\frac{p}{8}\right) \right);$$

In this case, the delay-difference operator is stable. Then, by applying Proposition 3.1 we obtain $p = 16.297...$ and $\tau = 0.6518...$, giving a critical solution at $s^* = \pm i3.125$ with multiplicity $m = 2$. To analyze this solution, we compute:

$$\left. \frac{\partial f}{\partial s} \right|_{(s^*, p^*, \tau^*)} = 0, \quad \left. \frac{\partial^2 f}{\partial s^2} \right|_{(s^*, p^*, \tau^*)} = 2.56516 + 3.63221i,$$

$$\left. \frac{\partial^\rho f}{\partial \tau^\rho} \right|_{(s^*, p^*, \tau^*)} = -0.638937 - 6.38911i.$$

Since $\frac{\partial f}{\partial \tau} \neq 0$ we have $\rho = 1$, and by Proposition 3.4-(i), there is only one ramification $s_{1,\sigma}$ (13) with $\beta_1 = \frac{1}{2}$ and $c_{1,\sigma}$, thus its asymptotic behavior around $s^*$ is described by:

$$s_{1,\sigma}(\tau) = s^* \pm (1.64 + 0.43i)(\tau - \tau^*)^{\frac{1}{2}} + o\left(|\tau - \tau^*|^{\frac{1}{2}}\right), \sigma = 1, 2.$$

Figure 3 illustrates the behavior of all solutions close to $s^*$.

## V. CONCLUDING REMARKS

In this paper, the asymptotic behavior of simple and multiple critical roots for quasi-polynomials of neutral type with commensurate delays has been considered. Specifically, we have proposed tools to analyze the asymptotic behavior of multiple critical solutions for quasi-polynomials of neutral type as a function of the system's parameters. The effectiveness of the proposed method is illustrated through some numerical examples. To the best of the authors' knowledge, such a characterization represents a novelty in the open literature.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Behavior of $s_{1,\sigma}$ around $s^* = \pm i3.125$.</figcaption></figure>

## ACKNOWLEDGMENTS

The work of C.-F. Méndez-Barrios has been supported by CONAHCyT-México, under the grant *Ciencia Básica y de Frontera: CBF2023-2024-2696*.

## REFERENCES

[1] G. Stépán, *Retarded Dynamical Systems: Stability and Characteristic Functions*. Longman Sci. and Tech., 1989.

[2] J. K. Hale and S. M. Verduyn Lunel, *Introduction to Functional Differential Equations*. New York: Springer, 1993.

[3] W. Michiels and S.-I. Niculescu, *Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach* (2. ed.), vol. 27 of Advances in design and control. Philadelphia: SIAM, 2014.

[4] V. Răsvan, "Dynamical systems with lossless propagation and neutral functional differential equations," in *Proc. MTNS (Padova, Italy*, pp. 527-532, 1998.

[5] V. Răsvan, "Propagation, delays, and stability (robust versus fragile)," *Int. J. Robust Nonlin. Contr.*, vol. 32, no. 6, pp. 3225-3250, 2022.

[6] A. Bellen, N. Guglielmi, and A. Ruehli, "Methods for linear systems of circuit delay differential equations of neutral type," *Circuits and Systems I: Fundamental Theory and Applications, IEEE Trans.*, vol. 46, pp. 212-215, 1999.

[7] K. Gopalsamy, X.-Z. He, and L. Wen, "On a periodic neutral logistic equation," *Glasg. Math. J.*, vol. 33, pp. 281-286, 1991.

[8] A. Martínez-González, C.-F. Méndez-Barrios, S.-I. Niculescu, J. Chen, and L. Félix, "Weierstrass approach to asymptotic behavior characterization of critical imaginary roots for retarded differential equations," *SIAM J. Control Optim.*, vol. 57, no. 1, pp. 1-22, 2019.

[9] P. Fu, S.-I. Niculescu, and J. Chen, "Stability of linear neutral time-delay systems: exact conditions via matrix pencil solutions," *IEEE Trans. Autom.*, vol. 51, no. 6, pp. 1063-1069, 2006.

[10] K. Gu, J. Chen, and V. L. Kharitonov, *Stability of time-delay systems*. New York: Springer Science & Business Media, 2003.

[11] N. Olgac and R. Sipahi, "A practical method for analyzing the stability of neutral type LTI-time delayed systems," *Automatica*, vol. 40, no. 5, pp. 847-853, 2004.

[12] C.-F. Méndez-Barrios, S.-I. Niculescu, and J. Chen, "Output feedback stabilisation of single-input single-output linear systems with i/o network-induced delays. an eigenvalue-based approach," *International Journal of Control*, vol. 87, no. 2, pp. 346-362, 2014.

[13] A. Martínez-González, S.-I. Niculescu, J. Chen, C. F. Méndez-Barrios, J. G. Romero, and G. Mejía-Rodríguez, *Delays and Interconnections: Methodology, Algorithms and Applications*, ch. Asymptotic Analysis of Multiple Characteristics Roots for Quasi-polynomials of Retarded-Type. Cham: Springer International Publishing, 2019.

[14] L. Berezansky and E. Braverman, "Asymptotic properties of neutral type linear systems," *Journal of Mathematical Analysis and Applications*, vol. 497, no. 2, p. 124893, 2021.

[15] C. T. C. Wall, *Singular points of plane curves*. New York: Cambridge University Press, 2004.

[16] J. Chen, G. Gu, and C. N. Nett, "A new method for computing delay margins for stability of linear delay systems," *Syst. Control Lett.*, vol. 26, p. 107-117, sep 1995.