The generic multiplicity-induced-dominancy property from retarded to neutral delay-differential equations: When delay-systems characteristics meet the zeros of Kummer functions

In this paper, which is a direct continuation and generalization of the recent works by the authors [https://doi.org/10.1051/cocv/2019073, https://doi.org/10.1016/j.jde.2021.03.003], we show the validity of the generic multiplicity-induced-dominancy property for a general class of linear functional differential equations with a single delay, including the retarded as well as the neutral cases. The result is based on an appropriate integral representation of the corresponding characteristic quasipolynomial functions involving some appropriate degenerate hypergeometric functions.

Given k, n ∈ N with k ≤ n, the binomial coefficient n k is defined as n k = n! k!(n−k)!and this notation is extended to k, n ∈ Z by setting n k = 0 when n < 0, k < 0, or k > n.For α ∈ C and k ∈ N, (α) k is the Pochhammer symbol for the ascending factorial, defined inductively as (α) 0 = 1 and (α) k+1 = (α + k)(α) k .
Finally, for the sake of simplicity in the formulations, we consider that the indices of rows and columns of matrices start from 0. More precisely, given n, m ∈ N * , an n × m matrix A is described by its coefficients A jk for integers j, k with 0 ≤ j < n and 0 ≤ k < m.

Introduction
This paper addresses the asymptotic behavior of the generic delay-differential equation (DDE): where the unknown function y is real-valued, n is a positive integer, m is a nonnegative integer such that m ≤ n, a k , α l ∈ R for k ∈ 0, n − 1 and l ∈ 0, m are constant coefficients, and τ > 0 is a delay.Under the assumption that α m = 0, equation (1.1) is a delay-differential equation which is said of retarded type if m < n and of neutral type if m = n (see, e.g., [28] and references therein).
For linear dynamical systems including delays in their model representation, spectral methods can be used to understand the asymptotic behavior of solutions by considering the roots of some characteristic function (see, e.g., [8,23,28,40,54,59]) which, for (1.1), is the function ∆ : C → C defined for s ∈ C by More precisely, the exponential behavior of solutions of (1.1) is given by the real number γ 0 = sup{ (s) | s ∈ C, ∆(s) = 0}, called the spectral abscissa of ∆, in the sense that, for every ε > 0, there exists C > 0 such that, for every solution y of (1.1), one has |y(t)| ≤ Ce Moreover, all solutions of (1.1) converge exponentially to 0 if and only if γ 0 < 0. An important difficulty in the analysis of the asymptotic behavior of (1.1) is that, contrarily to the delay-free case, the corresponding characteristic function ∆ has infinitely many roots.In all generality, the problem of characterizing the domain in the space of the equation's parameters that guarantee the exponential stability of solutions is a question of ongoing interest, see for instance [39].Since the 1950s, the complete understanding of the first-order retarded scalar differential equation with a single delay has benefited from the seminal work by Hayes [30], where the author gives a complete characterization of the rightmost spectral value location as a function of the delay and system's parameters.More recently, several works exploited Hayes results in control problems such as in delayed feedback and in stabilization problems.Unfortunately, Hayes approach remains complicated and natural extensions to nonscalar or higher-order retarded or neutral DDEs do not exist.
In control theory, the first pole-placement paradigm for time-delay systems, called finite pole placement, was introduced in the late 1970s in [34,44], where a prediction of the state over a delay interval is set to counteract the effect of the delay, hence reducing the closedloop system to a finite-dimensional plant.For an exhaustive presentation of the main ideas as well as some comparisons with other methods devoted to the control of dynamical systems including input delays, we refer to [57].Later, such an idea was deeply investigated in a more appropriate algebraic setting by [20] thanks to the introduction of the ring E defined as the set of all the meromorphic functions in the complex plane C that are of the form of P (s, e −h s )/Q(s), where Q is a polynomial in the Laplace complex variable s, P is a bivariate polynomial in s and e −h s , and h is a fixed positive real number.Indeed, the major issue for the algebraic design of controllers of linear time-invariant differential time-delay systems is the algorithmic study of the ring E. Furthermore, in practice, the limitation of such a paradigm was remarked during the early 2000s in [26] through the numerical design of a delayed controller theoretically able to stabilize a dynamical system described by a first-order scalar differential equation but leading to a closed-loop system whose stability is extremely sensitive to infinitesimal uncertainties.Such an instability mechanism, which is known as the spillover problem, is explained by the infinitely many spectral values generated by the discretization of the predictor.Furthermore, these fortuitous spectral values dominate in part the assigned spectral value.
Pole placement for delay systems goes beyond of a quasipolynomial interpolation problem.For instance, it has been shown in the 1980s by Ackermann [1] that n poles of the system can be assigned to (some) desired positions in the complex plane by n feedback parameters in the same way as in the finite-dimensional case.However, it is well known that such an interpolation is an effective placement if and only if the remaining spectral values of the closed-loop system are located to the left of the assigned poles, thus avoiding the spillover problem.In other words, such a strategy works if the assigned poles are dominant.But, as pointed out by [56], such a property is in general not guaranteed as commented in [51], where the proposed pole placement relies on rather heuristic trial-anderror placement of the dominant poles, making the pole placement procedure repeated with a different selection of assigned poles without any attempt to prove analytically the corresponding property.
To the best of our knowledge, the first "automated" pole placement for retarded timedelay systems is the numerical paradigm known as continuous pole placement introduced in [38].Unlike finite spectrum assignment, continuous pole placement does not render the closed-loop system finite-dimensional, but consists instead of controlling the corresponding rightmost eigenvalues.Such an idea represents a simple generalization of the pole placement for finite-dimensional systems represented by ordinary differential equations.It is based on the continuous dependence of the characteristic roots on the controller parameters, and the control strategy can be summarized as follows: "Shift" the unstable characteristic roots from C + to C − in a "quasi-continuous" way subject to the strong constraint that, during this shifting action, stable characteristic roots are not crossing the imaginary axis from C − to C + .We refer to [40] and references therein for further insights on the number of controlled characteristic roots (which is related to the available degrees of freedom induced by the controller structure) as well as the interpretation of continuous pole placement as a local strategy to solve an appropriate optimization problem where the objective function (rightmost root) is not differentiable.It is worth mentioning that continuous pole placement, initially applied to delay systems of retarded type, was extended to neutral systems in [41].
A more recent pole placement analytical paradigm ensues from an observation on the effect of multiple spectral values on the stability of DDEs, a property called multiplicityinduced-dominancy (MID) in [17].Indeed, some works have shown that, for some classes of time-delay systems, a real root of maximal multiplicity is necessarily the rightmost root, a property we call generic multiplicity-induced-dominancy, or GMID for short.This link between maximal multiplicity and dominance has been suggested in [46] after the study of some simple, low-order cases, but without any attempt to address the general case.To the best of the authors' knowledge, very few works have considered this question in more details until recently in works such as [9,[17][18][19][35][36][37]52].These works consider only DDEs with a single delay and show either the MID or the GMID property for each system under consideration.For instance, the MID or GMID properties are shown to hold for retarded equations of order 1 in [19], which proves dominance by introducing a factorization of ∆ in terms of an integral expression when it admits a root of maximal multiplicity 2; for retarded equations of order 2 with a delayed term of order zero in [18], using also the same factorization technique; or for retarded equations of order 2 with a delayed term of order 1 in [17], where both the MID and the GMID properties are investigated, using Cauchy's argument principle to prove dominance of the multiple root.Most of these results are actually particular cases of a more general result on the GMID property from [35] for generic retarded DDEs of order n with delayed "term" (polynomial) of order n − 1 (i.e.m = n − 1 in (1.1)), which relies on links between quasipolynomials with a real root of maximal multiplicity and the Kummer confluent hypergeometric function in terms of the location of the characteristic roots.The GMID property was also extended to neutral DDEs of orders 1 and 2 in [9,33,36], as well as to the case of complex conjugate roots of maximal multiplicity in [37].
The fact that a spectral value achieves maximum multiplicity imposes algebraic constraints on each of the system's "entries" (polynomial coefficients as well as the delay parameter).An MID-based approach is proposed in [4,5] operating the intimate representation of the quasipolynomial to provide conditions for one spectral value with an eligible intermediate multiplicity.This makes it possible to split the system parameters into two categories, some of them considered as model parameters (assumed to be fixed and known) and the remaining ones considered as values to be adjusted.Such a classification opens interesting perspectives in control design, such as the systematic tuning of the gains of the well-known Proportional-Integral-Derivative (PID) controller, able to stabilize singleinput/single-output plants including one delay in the input/output channel, as suggested in [33].
The contributions of this paper are fourfold.First, we propose an unified analytical framework for the characterization of the generic multiplicity-induced-dominancy of dynamical systems represented by DDEs in both retarded and neutral cases.The proposed methodology makes use of some degenerate hypergeometric functions and the existing links between such functions and the factorization of the characteristic functions of DDEs.In this sense, a first step was proposed by the same authors in [35], where only the special case of retarded equations with m = n − 1 was treated by using a result by Wynn [60] on the location of zeros of some Kummer confluent hypergeometric functions, the latter being proved using a representation of a quotient of Kummer functions as continued fractions.However, the underlying methodology from [60] seems hard to extend to more general retarded cases or to the neutral one, and for this reason the present paper is based on a different approach.
Second, our method is constructive and allows a better understanding of the existing links between the location of the characteristic root with generic maximal multiplicity, the so-called rightmost root, i.e., the spectral abscissa of the corresponding system, and the parameters of the dynamical systems.In particular, in some cases, it may be useful to minimize the spectral abscissa with respect to some of the tuning parameters with a guarantee of the exponential stability of the zero solution if the characteristic root satisfying the GMID property is located in C − .
Third, in terms of control, as mentioned above, the paradigm of generic multiplicityinduced-dominancy opens the perspective to a new control method, the so-called partial pole placement.More precisely, the GMID idea may be particularly adapted for tuning low-complexity controllers, i.e., controllers including a small number of parameters (including also the delay among the parameters) with a guarantee on the location of the remaining characteristic roots for the closed-loop system.For illustrating such an idea, some examples of PID controllers including a delay in the input/output channel are proposed for low-order systems in both retarded and neutral cases.
Finally, the ideas above can be tested and illustrated by using the software "Partial pole placement via delay action", or P3δ1 for short, which offers the computation of the values of the free parameters (typically tuning the controller gains) ensuring a prescribed multiplicity, establishes a certified assignment region for the rightmost root, performs a numerical computation to illustrate the distribution of the quasipolynomial characteristics, offers the possibility of a numerical study of the sensitivity of the spectrum with respect to uncertain parameters variations, and is also able to perform appropriate time-domain simulations, see for instance [14].New functionalities are now available covering both retarded and neutral cases.The remaining of the paper is organized as follows.Section 2 presents some prerequisites in complex analysis.It starts by recalling some qualitative properties of quasipolynomials and states some technical results on the distribution of zeros of some degenerate hypergeometric functions.The main results are presented in Section 3, where the GMID property is proved for the general retarded and neutral differential equations with a single delay.These results, which are implemented in the control design software P3δ, helped us to study standard examples in the design of stabilizing proportional-integral-derivative controllers as in [33].Section 4 provides some illustrations of the ideas above applied to a couple of comprehensive control problems.More precisely, the stabilization of the classical pendulum and the feedback stabilization for a scalar conservation law with PI boundary control, a description of the new features of the P3δ software as well as some further remarks on the optimization of the trivial solution decay are presented.Some concluding remarks end the paper.

Preliminaries and prerequisites
This section recalls some preliminary results on properties of quasipolynomial functions (Section 2.1), the distribution of zeros of Laplace transform of Bernstein polynomials (Section 2.2), and properties and the distribution of nonasymptotic zeros of degenerate hypergeometric functions (Section 2.3).

Quasipolynomials
Let us start by recalling the classical definition of a quasipolynomial, also known as exponential polynomial, and its degree (see, e.g., [10,58]).
A quasipolynomial Q is an entire function Q : C → C which can be written under the form where is a nonnegative integer, σ 0 , . . ., σ are pairwise distinct real numbers, and, for k ∈ 0, , P k is a nonzero polynomial of degree d k ≥ 0 with coefficients belonging to C. The integer D P S = + k=0 d k is called the degree of Q.When σ 0 = 0 and σ k < 0 for k ∈ 1, in the above definition, Q is the characteristic function of a linear time-delay system with delays −σ 1 , . . ., −σ , which motivates our study of such quasipolynomials.Notice that, if the coefficients of a quasipolynomial Q are real, which is the case in (1.2), then the corresponding zeros are necessarily either real or complex-conjugate pairs.The roots of a quasipolynomial do not change when its coefficients are all multiplied by the same nonzero number, and hence one may always assume, without loss of generality, that one nonzero coefficient of a quasipolynomial is normalized to 1, such as the coefficient of the term of highest degree in P 0 , which is the case in (1.2), for instance.The degree of a quasipolynomial is the number of the remaining coefficients.Since each polynomial P k of degree d k has d k + 1 coefficients, the quasipolynomial degree is then the sum of these numbers discounting the normalized coefficient, giving rise to the number It is worth to note that a quasipolynomial admits an infinite number of roots, except in trivial cases when the quasipolynomial reduces to a polynomial.Fortunately, there does exist a link between the degree of a quasipolynomial and the number of its roots in horizontal strips of the complex plane, thanks to a classical result known as Pólya-Szegő bound (see [48, Part Three, Problem 206.2]), which we state in the next proposition.
Proposition 2.1.Let Q be a quasipolynomial of degree D P S given under the form (2.1), α, β ∈ R be such that α ≤ β, and Notice that the above result was first introduced and claimed in the problems collection published in 1925 by G. Pólya and G. Szegő.In the fourth edition of their book [48, Part Three, Problem 206.2],G. Pólya and G. Szegő emphasized that the proof was obtained in the meantime by N. Obreschkoff using the principle argument, see [43].
As an immediate consequence of Proposition 2.1, given a root s 0 ∈ C of a quasipolynomial Q of degree D P S , by letting β = α = (s 0 ) in the statement of Proposition 2.1, one concludes that s 0 has multiplicity at most D P S .Hence, for a quasipolynomial ∆ under the form (1.2), which is of degree m + n + 1, any of its roots has multiplicity at most m + n + 1.
Since quasipolynomial functions admits an infinite number of zeros, an important information on the distribution of such zeros is the location of the corresponding dominant or rightmost root.Let us consider a function Q : C → C and a value s 0 ∈ C such that Q(s 0 ) = 0. We say that s 0 is a dominant (respectively, strictly dominant) root of Q if, for every s ∈ C \ {s 0 } such that Q(s) = 0, one has (s) ≤ (s 0 ) (respectively, (s) < (s 0 )).
The next lemma, which is a key ingredient to simplify the proofs of our main results in Section 3, describes how the coefficients of a quasipolynomial change under a translation and a dilation of the complex plane.Its proof, omitted here, can be carried out by straightforward computations and is very similar to that of [35,Lemma 4.1].Lemma 2.2.Let s 0 ∈ R, ∆ be the quasipolynomial from (1.2), and consider the quasipolynomial ∆ : C → C obtained from ∆ by the change of variables z = τ (s − s 0 ) and multiplication by τ n , i.e., where (2.4) The relations between the coefficients b 0 , . . ., b n−1 , β 0 , . . ., β m and a 0 , . . ., a n−1 , α 0 , . . ., α m expressed in (2.4) can be expressed under matrix form as .

Zeros' distribution of finite Laplace transform
We discuss in this section Laplace transforms of some integrable functions f concentrated on a finite support, which, without loss of generality, can be taken as the interval (0, 1), i.e., Before turning to the core of this section, we present the following result on the integral of the product of a polynomial and an exponential, which is rather simple but of crucial importance to prove our main result.Its proof is straightforward and can also be found in [ The question of locating the zeros of functions F under the form (2.6) is not new and was for instance the subject of several pioneering works by Hardy [29], Pólya [47], and Titchmarsh [55].Furthermore, beyond its strict mathematical importance, the location of zeros of (2.6) is related to a wide range of problems coming from physics and engineering.Despite having been established more than a century ago, the following result from [47] remains relevant and can be found in [49, Part Five, Chapter 3, Problem 177, page 66] (see also [53] for further related results).The assumptions on f in Theorem 2.5 (namely, positivity and monotonicity of the function f ) are quite restrictive for our purpose.Similar results in locating the zeros of (2.6) can be obtained when f is a Bernstein polynomial [11], that is, see also [6] for further insights on Bernstein polynomials.The following simple example provides an illustration of that fact.
Example 2.6 (Zeros' frequency bound approach).Consider the function F : C → C defined for s ∈ C by We will show that all zeros of F lie in the open right half-plane C + .Notice first that, using Proposition 2.4, one obtains Obviously, F is an analytic function since the numerator of the right-hand side of (2.9) admits a root at zero with multiplicity four.Noticing that this numerator is a quasipolynomial of degree 4, we deduce, as a consequence of Proposition 2.1, that the function F does not admit any real root.One can also show that no zero of F lies on the imaginary axis.As a matter of fact, assume that there exists a root s 0 of F such that s 0 = i ω 0 with ω 0 = 0. Then The vanishing of the corresponding real and imaginary parts allows to eliminate the trigonometric functions as rational expressions of the crossing frequency ω 0 .Namely, one obtains Further, using the standard property sin 2 (ω 0 ) + cos 2 (ω 0 ) = 1, one deduces that ω 0 must satisfy ω 4 0 (ω 2 0 + 9) = 0, which is impossible since ω 0 = 0, yielding the required conclusion.We are left to show that F has no zeros with negative real part.To investigate the potential roots of F with negative real parts we rather investigate the roots of s → F (−s) with positive real parts.Indeed, define the auxiliary function (2.10) Obviously, any nontrivial zero s 0 = ξ + i ω of G satisfies F (−s 0 ) = 0. We will first show that any such s 0 with positive real part must necessarily satisfy 0 < |ω| < π.
Since ξ is assumed to be positive, we have e 2ξ > 1 + 2ξ and then one deduces from equation (2.10) that the pair (ξ, ω) necessarily satisfies which, using the notation Ω = ω 2 , is equivalent to H ξ (Ω) < 0, where Since We can now use the fact that any root s 0 = ξ + iω of G with ξ > 0 satisfies 0 < |ω| < π to deduce that no such root can exist.Indeed, taking the imaginary part of the equality F (−s 0 ) = 0 and using (2.8), we get 0 = The fact that |ω| < π implies that the integrand is positive, yielding the result.Hence all roots s of (2.9) satisfy (s) > 0.
The location of roots of functions of type (2.7) with indices m and n which are not necessarily integers remains possible thanks to the existing link between Bernstein polynomials with degenerate hypergeometric functions, which is the topic of our next session.

Degenerate hypergeometric functions
As it will be proved in Section 3.1, when the quasipolynomial ∆ from (1.2) admits a root of maximal multiplicity m + n + 1, it can be represented in terms of a Kummer confluent hypergeometric function, whose nontrivial roots coincide with the nontrivial roots of another special function, known as Whittaker function.These families of special functions have been extensively studied in the literature, with, in particular, a wide range of results on the properties of the asymptotic distribution of their zeros (see, e.g., [21], [27, Chapter VI], [45,Chapter 13]).For the purposes of our paper, however, we need information on the location of all zeros of such functions, and not only their asymptotic distribution, a subject that has been considered in fewer works in the literature (we refer the interested reader to the introduction of [13] for a more detailed description on these questions).This section provides a brief presentation of the results that shall be of use in the sequel.We start by recalling the definition of Kummer confluent hypergeometric functions used in this paper.
where we recall that, for α ∈ C and k ∈ N, (α) k is the Pochhammer symbol (defined in the "Notation" section).which has a regular singular point at z = 0 and an irregular singular point at z = ∞.It is well known that (2.12) admits two linearly independent solutions, both of them usually called Kummer confluent hypergeometric functions.For our purpose we are concerned only with the solution (2.11).
We shall need the following classical integral representation of Φ, which can be found, for instance, in [21,25,27,45].Proposition 2.9.Let a, b ∈ C and assume that (b) > (a) > 0.Then, for every z ∈ C, we have where Γ denotes the Gamma function.
We also recall (see, e.g., [45, (13.2.39)]), that, for every a, b, z ∈ C such that b is not a nonpositive integer, we have The following result, which is proved in [13] using the Green-Hille transformation from [31], gives insights on the distribution of the nonasymptotic zeros of Kummer hypergeometric functions with real arguments a and b.The techniques used in that reference do not rely on Hille's approach, but use instead a continued fraction representation of a ratio of Kummer functions (see [60, Theorem B(ii)]).Note that Proposition 2.10 does not cover all cases of Wynn's result, since the assumption b ≥ 2 is equivalent to α ≥ 1 2 , whereas Wynn's result is obtained for all α > − 1 2 .

Main results
This section is dedicated to our main results on the location of zeros of the characteristic function ∆ given in (1.2) of system (1.1).Hereafter, we will characterize the manifold corresponding to the existence of a spectral value of (1.2) with multiplicity reaching the Pólya-Szegő bound D P S = m + n + 1 and will prove that such spectral value is necessarily dominant.

Maximal multiplicity and quasipolynomial factorization
The main result of this section is the following characterization of the existence of a real root attaining the maximal multiplicity D P S .
Theorem 3.1.Consider the quasipolynomial ∆ given by (1.2) and let s 0 ∈ R. The number s 0 is a root of multiplicity Remark 3.3.Under conditions (3.1), s 0 is the unique real root of ∆, and, more precisely, it is the unique root of ∆ on the horizontal strip {s ∈ C | | (s)| < 2π τ } of the complex plane.Indeed, applying Proposition 2.1 with α = 0 and any β ∈ (0, 2π τ ), we deduce that the number of roots of ∆ in the strip S β = {s ∈ C | 0 ≤ (s) ≤ β}, counted according to their multiplicity, is at most D P S , since τ (β−α) 2π < 1 in this case.Since s 0 ∈ S β and its multiplicity as a root of ∆ is exactly D P S , this shows that no other root of ∆ can belong to S β .The conclusion now follows since β ∈ (0, 2π τ ) is arbitrary, and by remarking that ∆ has only real coefficients and thus its roots appear in complex-conjugate pairs.Remark 3.4.Theorem 3.1 deals only with a real root s 0 of ∆ attaining the maximal multiplicity D P S .The main reason for this restriction is that, as shown in [16, Corollaries 1 and 2], nonreal roots of ∆ cannot have a multiplicity equal to the Pólya-Szegő bound D P S (see also [37] for further discussion on the maximal multiplicity of complex roots in the case n = 2 and m = 1).Indeed, any root s 0 of ∆ attaining the maximal multiplicity D P S necessarily satisfies (3.2), and thus it will be real since a n−1 is real.
The proof of Theorem 3.1 can be simplified by noticing that s 0 is a root of the quasipolynomial ∆ of a given multiplicity if and only if 0 is a root with the same multiplicity of the quasipolynomial ∆ defined by ∆(z) = τ n ∆(s 0 + z τ ).This transformation of ∆ into ∆ corresponds to the linear change of variable z = τ (s − s 0 ), which shifts the desired multiple root s 0 to 0 and reduces the delay τ to 1, and it has been described in Lemmas 2.2 and 2.3 in Section 2.1.We can then focus on providing necessary and sufficient conditions on the coefficients of ∆ in order for 0 to be a root of maximal multiplicity D P S = m + n + 1, which is the main subject of our next result, Lemma 3.5, which also states that, under such conditions, ∆ can be factorized as the product of z m+n+1 and an entire function expressed as the Laplace transform of a Bernstein polynomial, as the one discussed in Example 2.6.
Proof.We first prove that the right-hand side of (3.4) is indeed a quasipolynomial of the form (2.3) with coefficients b 0 , . . ., b n−1 , β 0 , . . ., β m given by (3.3) and admitting 0 as a root of multiplicity D P S .For that purpose, let us introduce the function Q : C → C defined as the right-hand side of (3.4), i.e., for every z ∈ C, Clearly, Q is an entire function and 0 is a root of multiplicity m+n+1 of Q.Let p : R → R be the Bernstein polynomial defined for t ∈ R by p(t) = t m (1 − t) n , whose degree is m + n.
In order to conclude the proof, it suffices to show that (3.3) is the unique choice of coefficients for ∆ ensuring that 0 is a root of multiplicity m + n + 1.To see that, notice that, since the degree of the quasipolynomial ∆ is m + n + 1, 0 is a root of multiplicity m + n + 1 of ∆ if and only if ∆ (k) (0) = 0 for every k ∈ 0, m + n .By [15,Proposition 5.1], this is the case if and only if the coefficients b 0 , . . ., b n−1 , β 0 , . . ., β m satisfy The first two equalities in (3.6) allow us to compute b 0 , . . ., b n−1 once β 0 , . . ., β m are known, and the last two equalities in (3.6) provide a linear system on β 0 , . . ., β m .This linear system can be written under the form T β = −e 0 , where T = (T j,k ) j,k∈ 0,m is defined by , and e 0 = (1, 0, . . ., 0) T ∈ R m+1 .By the previous arguments, the coefficients b 0 , . . ., b n−1 , β 0 , . . ., β m given by (3.3) necessarily satisfy (3.6), since they ensure that 0 is a root of multiplicity m + n + 1 of ∆, and thus we are left to prove that (3.6) admits a unique solution, which is equivalent to showing that the matrix T is invertible.Let us introduce the diagonal matrices D 0 = diag(0!, 1!, . . ., m!) and D n = diag(n!, (n+ 1)!, . . ., (m + n)!) and set B = D n T D −1 0 .Writing B = (B j,k ) j,k∈ 0,m , we have We have B = LU , where L = (L j,k ) j,k∈ 0,m and U = (U j,k ) j,k∈ 0,m are defined, for j, k ∈ 0, m , by where we use the identity m =0 j n k− = n+j k (see, e.g., [35,Proposition 2.16]).Notice that L is lower triangular and U is upper triangular, and hence the factorization B = LU corresponds to the LU factorization of B. As a consequence of this factorization, we also deduce that det B = (−1) n (m+1) and, in particular, B is invertible.Hence, T = D −1 n BD 0 is also invertible, concluding the proof.
Proof of Theorem 3.1.Define ∆ from ∆ as in (2.2).One immediately verifies that s 0 is a root of multiplicity m + n + 1 of ∆ if and only if 0 is a root of multiplicity m + n + 1 of ∆.The result then follows by a straightforward computation using Lemmas 2.3 and 3.5.

Dominance of roots of maximal multiplicity
Theorem 3.1 provides necessary and sufficient conditions for a real number s 0 to be a root of maximal multiplicity of the quasipolynomial ∆ from (1.2).The main result of this section states that, under those conditions, s 0 is necessarily a dominant root of ∆.
To prove Theorem 3.6, we rely, as in the proof of Theorem 3.1, on the linear change of variable z = τ (s − s 0 ).We thus first study the quasipolynomial ∆ from (2.3) under conditions (3.3), in which case we have, by Theorem 3.1, the factorization (3.4).The factorization (3.4) can also be written, thanks to (2.13), as where Φ is the Kummer confluent hypergeometric function defined in (2.11).The next lemma uses properties of the roots of Φ in order to deduce that 0 is a dominant root of ∆.
In order to obtain the additional characterization of the spectrum in the neutral case stated in Theorem 3.6(b), we also provide the following result on the location of the roots of ∆ in the case m = n.Lemma 3.9.Let n ∈ N * , m = n, b 0 , . . ., b n−1 , β 0 , . . ., β n ∈ R be given by (3.3), ∆ be the quasipolynomial given by (2.3), Ξ n be as in the statement of Theorem 3.6(b), and Proof.From (2.3) and (3.3), we have Hence, z = iζ is a root of ∆ if and only if The right-hand side of the above equality is equal to the complex conjugate of the left-hand side.Hence, the above equality is equivalent to which happens if and only if This last equality is equivalent to ζ ∈ Ξ n , concluding the proof.
Using Lemmas 3.8 and 3.9, we can easily conclude the proof of Theorem 3.6.
Proof of Theorem 3.6.Define ∆ from ∆ as in (2.2).For n > m, item (a) can be shown by noticing that, if s is a root of ∆ with s = s 0 , then, by (2.2), z = τ (s − s 0 ) is a root of ∆ with z = 0. Hence, by Lemma 3.8, (τ (s − s 0 )) < 0, showing, since τ > 0, that (s) < s 0 .Similarly, for m = n, the first part of item (b) can be shown by applying first Lemma 3.8, which gives (s) = s 0 for any root s of ∆, and the last part of item (b) follows by applying Lemma 3.9.

Consequences on stability
The third main result we present in this paper is the following one on the the stability of the trivial solution of (1.1), which is an immediate consequence of (3.2), Theorem 3.1, and Theorem 3.6.

Illustrative examples
PID controllers have been extensively used to control and regulate industrial processes which are typically modeled by reduced-order dynamics.In this section, we shall illustrate how one can tune the controller gains using the GMID property through a couple of comprehensive examples.The first case deals with a retarded equation, while the second one corresponds to a neutral equation.

Stabilizing the classical pendulum via delay action
Consider the dynamical system modeling a friction-free classical pendulum.The adopted model is studied in [3] and, in the sequel, we keep the same notations.The dynamics of a controlled classical pendulum are governed by the following second-order differential equation θ(t) + g L sin(θ(t)) = u(t) ( where θ(t) stands for the angular displacement of the pendulum at time t with respect to the stable equilibrium position, L is the length of the pendulum, g is the gravitational acceleration, and u(t) represents the control law, which stems from an applied external torque.The problem we consider is to make the open-loop stable equilibrium locally asymptotically stable via the action of the external torque u.This suggests first to linearize (4.1) around the trivial equilibrium and get Next assume that such a system is controlled using a standard delayed PD controller of the form The local stability of the closed-loop system is reduced to study the location of the spectrum of the quasipolynomial which is a quasipolynomial of degree 4. Remarks 3.3-3.4allow concluding that the maximal multiplicity 4 can be reached only by a root of ∆ on the real axis.
In this case, we apply the GMID property to establish the gains of the stabilizing delayed PD controller.As a matter of fact, in this case, the only admissible quadruple root is given by which is achieved if the controllers' gains and the delay τ are taken as With the proposed choice of the gains and the delay, conditions (3.1) are satisfied for the quasipolynomial (4.4), and thus, by Theorem 3.6, s 0 given by (4.5) is the spectral abscissa of the closed-loop system (4.2)-(4.3),ensuring then the exponential stability of the trivial solution as announced in Theorem 3.10.Remark 4.1.Notice that the GMID consists in forcing a root to reach its maximal multiplicity, which does not allow any degree of freedom in assigning s 0 , as precised in (4.5).In order to allow for some additional freedom when assigning s 0 , one can relax such a constraint by forcing the root s 0 to have a multiplicity lower than the maximal, and also consider the delay as a free tuning parameter.This motivates the study of a (non-generic) MID property, which was carried out, for instance, in [17] for second-order systems.For instance, the only admissible triple roots of (4.4) are given by which are achieved if τ < 2 L/g and the controllers' gains are taken as Following [17, Theorem 4.2], the triple root at s = s + is the rightmost root of (4.4).Thus the delay, if seen as a tuning parameter, allows to assign the rightmost root at s = s + arbitrarily large (in absolute value) for small delay, as illustrated in Figure 4.1.

Feedback stabilization for a scalar conservation law with PI boundary control
As an illustrative example of the GMID for neutral equations, we consider the problem of stabilization of solutions of a partial differential equation of hyperbolic type.More precisely, we revisit the problem of exponential stabilization of the following scalar conservation law proposed in [24]: where L > 0 and ϕ(t, x) denotes the system state at position x ∈ (0, L) and in time t ∈ [0, +∞).As considered in [24], the value λ, which denotes the velocity of propagation, is assumed to be a positive constant.Equation (4.6) is accompanied with a boundary condition under the form of a PI controller: where k p and k i are the feedback parameters representing proportional and integral control gains.Applying the Laplace transform to both sides of the boundary condition and multiplying by s one obtains the closed-loop characteristic function which corresponds to the characteristic function of a first-order neutral equation, that is, a function under the form (1.2) with m = n = 1.In this case, the degree D P S of ∆ is equal to 3 and, as mentioned in Remarks 3.3-3.4,the maximal multiplicity can be achieved only by a root on the real axis.Next, by exploiting the results of Theorem 3.1, Theorem 3.6, and Theorem 3.10, we conclude that forcing a triple spectral value guarantees its dominance as a root of (4.8), and then the exponential stability of solutions of (4.6)-(4.7).More precisely, by tuning the controller gains as one achieves the unique admissible triple root, which is s 0 = − 2 λ L and corresponds to the decay rate of solutions of (4.6)-(4.7).Furthermore, as shown in Theorem 3.6(b), the set of roots of ∆ is s 0 shows the result of a numerical computation of the roots of (4.8) with the parameters (4.9), while

New features of P3δ
P3δ, whose name stands for Partial pole placement via delay action, is a Python software with a friendly user interface for the design of parametric stabilizing feedback laws with time delays, thanks to properties of the distribution of quasipolynomials' zeros.It exploits mainly the MID/GMID property as well as the coexisting real roots-induced-dominancy or CRRID for short, established in [2,7,14].Initially, P3δ has been established for the control design for retarded differential equations.Based on the results of this work as well as the references [9,33,36], the software has been updated and is now able to treat linear neutral functional differential equations as well.P3δ is freely available for download on https://cutt.ly/p3delta,where installation instructions, video demonstrations, and the user guide are also available.Interested readers may also contact directly any of the authors of the paper.Since its creation, P3δ had vocation to be available to the greatest number and on all possible platforms.The current version of the software is available in local executable version as well as an online version ready to use in one click.The online version of P3δ is hosted on servers thanks to the Binder service [50].Binder allows to create instances of personalized computing environment directly from a GitHub repository that can be employed and shared by users.The Binder service is free to use and is powered by BinderHub, an open-source tool that deploys the service in the cloud.The online version of P3δ is written in Python and structured as a Jupyter Notebook, an open document format which can contain live code, equations, visualizations, and text.The Jupyter Notebook is completed by a friendly user interface built using interactive widgets from Python's ipywidgets module.

Further remarks on the maximal damping
Although there exists a root locus approach for single delay time-delay systems [32] able to monitor the variation of the spectrum with respect to parameters variation, then allowing to characterize the spectral abscissa, unfortunately this method doesn't recover the multiple root case due to the lack of regularity with respect to parameter's variation.
To the best of the authors' knowledge, the problem of the selection of the system's parameters guaranteeing the maximum damping of reduced-order time-delay systems solutions goes back to the almost forgotten contributions of Pinney in the 1960s [46], where a complex analysis method based on Cauchy index theorem coupled with the D-partition method [42] are used.Furthermore, the link between a spectral value of maximal multiplicity and the maximum damping has been emphasized for first-and second-order DDEs, however without providing explicit proofs.As a matter of fact, it has been observed that the minimization of the spectral abscissa suggests the "exhaustion" of the equation parameters, which is achieved in Equation (3.1) from Theorem 3.1.
Next, in the context the design of output feedback able to achieve a fast stabilization of finite-dimensional systems, an optimization technique called Convex hull technique has been proposed in [22], where the existing link between the maximum damping and the characteristic root of maximum multiplicity has been established, see also [12].More recently, to recover such a link for second-order delay systems, an optimization approach has been used in [52].Notice that the extension of such optimization methods to neutral DDEs or even to retarded DDEs of fixed but higher orders (greater than 3) systems is a quite complicated question especially in the presence of free parameters.

Concluding remarks
In this paper, we have further investigated the multiplicity-induced-dominancy property for single-delay linear functional differential equations.Thanks to the reduction of the corresponding characteristic function to an appropriate integral representation that corresponds to some degenerate Kummer hypergeometric function, we have shown the validity of the generic multiplicity-induced-dominancy (GMID) property; that is, the characteristic spectral values of maximal multiplicity are necessarily dominant for retarded as well as neutral delay-differential equations of arbitrary order.
Notice that the MID property may hold even when it is about a spectral value with a strictly-intermediate admissible multiplicity, as shown in [17] for second-order plants and in [4] for n th -order retarded equations admitting a real spectral value with multiplicity n + 1 and a finite dimensional part admitting exclusively real modes.However, in all generality, the limits of the MID property remain an open question.
All these arguments make our constructive method based on confluent hypergeometric functions applied to linear delay-differential equations (retarded as well as neutral) of arbitrary order particularly attractive and advantageous in many ways and, in particular, in control design problems involving free parameters to be tuned.As shown through some illustrative examples concerning classical PID control schemes, the partial pole placement method that is based on the use of GMID property may appear as a good alternative to the existing control methodologies.Indeed, the corresponding controllers are of "lowcomplexity" (i.e., reduced number of parameters to be tuned) and there exist explicit guarantees on the location of the remaining characteristic roots if the GMID is reached.

Notation.
In this paper, N * denotes the set of positive integers and N = N * ∪ {0}.The set of all integers is denoted by Z and, for a, b ∈ R, we denote a, b = [a, b] ∩ Z, with the convention that [a, b] = ∅ if a > b.For a complex number s, (s) and (s) denote its real and imaginary parts, respectively.The open left and right complex halfplanes are the sets C − and C + , respectively, defined by C − = {s ∈ C | (s) < 0} and C + = {s ∈ C | (s) > 0}.

Theorem 2 . 5
(G.Pólya, 1918).Let f be a positive and continuously differentiable function defined in the interval [0, 1] and satisfying f (t) < 0 for every t ∈ [0, 1].Consider the function F : C → C defined from f as in(2.6).Then all the zeros of F lie in the open right half-plane C + .

Definition 2 . 7 .
Let a, b ∈ C and assume that b is not a nonpositive integer.The Kummer confluent hypergeometric function Φ(a, b, •) : C → C is the entire function defined for z ∈ C by the series

(3. 1 )Remark 3 . 2 .
Before turning to the proof of Theorem 3.1, a few remarks are in order.By considering the first equation in (3.1) with k = n − 1, one obtains the simple and interesting relation between s 0 , τ , and a n−1 given by

Theorem 3 . 6 .
Consider the quasipolynomial ∆ given by (1.2) and let s 0 ∈ R. (a) (Retarded) If m < n and (3.1) is satisfied, then s 0 is a strictly dominant root of ∆.(b) (Neutral) If m = n and (3.1) is satisfied then, s 0 is a dominant root of ∆ and, for every other complex root s of ∆, one has (s) = s 0 .More precisely, the set of roots of ∆ is {s 0

Figure 4 . 1 :
Figure 4.1: The behavior of the triple root (spectral abscissa) of (4.4) at s = s + as a function of the free delay parameter 0 < τ < 2 L/g for L/g = 1.
35, Proposition 2.1].Let d ∈ N and p be a polynomial of degree at most d.Then, for every z