# The Distance to the S147 Supernova Remnant

C. S. Kochanek Department of Astronomy, The Ohio State University
140 West 18th Avenue, Columbus, OH 43210  Center for Cosmology and
Astroparticle Physics, The Ohio State University
191 W. Woodruff Avenue, Columbus, OH 43210  John C. Raymond  Harvard-
Smithsonian Center for Astrophysics
60 Garden Street, Cambridge, MA 02138  Nelson Caldwell  Harvard-Smithsonian
Center for Astrophysics
60 Garden Street, Cambridge, MA 02138

###### Abstract

In the absence of a parallax distance to a pulsar or a surviving binary in a
supernova remnant (SNR), distances to Galactic SNRs are generally very
uncertain. However, by combining Gaia data with wide field, multi-fiber
echelle spectroscopy, it is now possible to obtain accurate distances to many
SNRs with limited extinction by searching for the appearance of high velocity
Ca II or Na I absorption lines in hot stars as a function of distance. We
demonstrate this for the SNR S147 using the spectra of 259 luminous, blue
stars. We obtain a median distance of $1.37$ kpc ($1.30$ to $1.47$ kpc at 90%
confidence) that is consistent with the median parallax distance to the pulsar
of $1.46$ kpc ($1.12$ to $2.10$ kpc at 90% confidence), but with significantly
smaller uncertainties. Our distance is also consistent with the distance to
the candidate unbound binary companion in this SNR, HD 37424. The presence of
high velocity absorption lines is correlated with the emission line flux of
the SNR but not with the radio flux.

Supernova remnants (1667) – interstellar line absorption (843) – distance
measure (395)

††facilities: MMT

##  1 Introduction

Supernova remnants (SNR) dominate the dynamics of the ISM and the chemical
evolution of the Galaxy. They are physically interesting and important probes
of both the explosion and the nature of the progenitor stars. For example,
estimates of ejecta masses and velocities constrain the progenitor and
explosions (e.g., Temim et al. 2022). They are laboratories for shock physics
(e.g., Raymond et al. 2020a, Raymond et al. 2020b). They can be searched for
bound or unbound binaries after the explosion in both core collapse (e.g.,
Kochanek et al. 2019, Kochanek 2021) and Type Ia remnants (e.g., Kerzendorf et
al. 2009). Just as in extragalactic SNRs (e.g., Badenes et al. 2009, Jennings
et al. 2014), the age distribution of the surrounding stellar population can
be used to estimate the mass of the progenitor (Kochanek 2022).

A major challenge in using SNRs for any of these applications is determining
their distances. Almost all Galactic SNR distances are significantly
uncertain. The most common estimates use HI velocities or pulsar dispersion
measures, neither of which are very precise. For our eventual target, S147
(Simeis 147, G180.0$-$1.7, the Spaghetti nebula), the traditional distance
estimates range from roughly 0.8 to 2 kpc (see the summary in Dinçel et al.
2015), leading to factor of 2.5 uncertainties in sizes or velocities from
proper motions and a factor of 6 in luminosities or masses. Only three
Galactic SNRs have distances from pulsar parallaxes: Vela (Dodson et al.
2003), the Crab (Lin et al. 2023, both VLBI and Gaia) and S147. The pulsar in
S147, PSR J0538$+$2817, has a parallax of $0.72\pm 0.12$ mas (Ng et al. 2007,
Chatterjee et al. 2009) which corresponds to a median distance of $1.46$ kpc
($1.12$-$2.10$ kpc at 90% confidence). There are also parallax distances to
G039.7$-$02.0, G205.5$+$00.5 (Monoceros loop) and G284.3$-$01.8 because they
contain the interacting binaries SS 433 (see the review by Margon 1984), MWC
148 (Hinton et al. 2009) and 2MASS J10185560$-$5856459 (Corbet et al. 2011),
respectively, and the donor stars have Gaia (Gaia Collaboration et al. 2016,
Gaia Collaboration et al. 2023) parallaxes. Curiously, the Gaia parallax for
SS 433 of $\varpi=0.118\pm 0.023$ mas (geometric distance $6.4$-$8.5$ kpc,
photogeometric distance of $5.5$-$8.9$ kc, Bailer-Jones et al. 2021) is at
best marginally consistent with estimates from kinematic models of its jets
(e.g., $5.5\pm 0.2$ kpc, Blundell & Bowler 2004).

Cha et al. (1999) demonstrated a new method for measuring SNR distances,
looking for the appearance of high velocity (up to several 100 km/s) Ca II and
Na I absorption features from the SNR in the spectra of hot stars (which lack
such features) superposed on the SNR as a function of distance. They estimated
a distance to the Vela SNR of $250\pm 30$ pc which was later confirmed by the
pulsar parallax distance of $286\pm 16$ pc (Dodson et al. 2003).
Interestingly, the high velocity absorption features are also rapidly time
variable in strength and wavelength (e.g., Kameswara Rao et al. 2020 most
recently). The only other case where this approach has been used is for the
Cygnus Loop (Fesen et al. 2018, Fesen et al. 2021, Ritchey et al. 2024).

Cha et al. (1999) depended on Hipparcos (Perryman et al. 1997) parallaxes and
so could only do a very nearby SNR like Vela. With Gaia, it is possible to
apply the method to tens of SNRs, as demonstrated by the Cygnus Loop result
(Fesen et al. 2018, Fesen et al. 2021). But both studies depended on observing
stars one at a time, which leads to either a very expensive program ($\sim 70$
stars for Vela) or a very sparse sample ($6$ in each study for Cygnus).
Interestingly, only some background stars show the high velocity absorption
features. It is not known if the existence of absorption is (anti)correlated
with any observational property of the SNR (radio, emission lines, X-rays,
etc.). Fortunately, there are now wide field, multiobject, fiber fed, echelle
spectrographs like Hectochelle on the MMT (Fabricant et al. 2005, Szentgyorgyi
et al. 2011), M2FS on Magellan (Mateo et al. 2012), FLAMES on the VLT
(Pasquini et al. 2000), and HERMES on the AAT (Sheinis 2016) which allow
observations of tens to hundreds of stars simultaneously.

Here we demonstrate this approach using the MMT and Hectochelle to observe
stars in the direction of S147. The one degree field of view of the fiber
positioner is well matched to the $3^{\circ}$ diameter (Green 2009) of S147
and many other SNRs. Hectochelle can obtain 240 spectra at a resolution of
34,000 ($\sim 8$ km/s) in 150Å windows centered on Ca II $\lambda\lambda
3933.663,3968.47$ (H and K) or Na I $\lambda\lambda 5889.51,5895.924$ (D1,
D2). We chose S147 as a demonstration case because the pulsar parallax (Ng et
al. 2007, Chatterjee et al. 2009) provides a check of the method. S147 is also
of interest because it contains the one fairly convincing example of a binary
star, HD37424, unbound in a supernova explosion (Dinçel et al. 2015, Kochanek
2021). We describe the selection of the stars and the observations in §2, and
the observational results and analysis in §3. In §4 we summarize the results
and discuss future applications.

Figure 1: (Left) The distribution of the stars projected on S147 in absolute
magnitude $M_{G}$ and $B_{P}-R_{P}$ color after correcting for extinction. The
vertical scale on the right is the corresponding apparent magnitude at the
parallax distance of the pulsar. The curves are PARSEC isochrones from
$10^{7}$ to $10^{9}$ years in steps of $0.5$ dex colored solid red (dashed
magenta) for temperatures $>10^{4}$ K ($<10^{4}$ K). The red square is the
candidate unbound binary companion, HD 37424, and the red triangles are two
background stars, HD 36665 and HD 37318, with known high velocity absorption
lines. Stars in the upper left corner outlined by the red box are our high
priority targets. The yellow/red pentagons are the two stars whose spectra are
shown in Figs. 6 (bluer star) and 6 (redder star).  Figure 2: (Right)
Heliocentric velocities of Ca II (black squares) and Na I (red triangles)
absorption features as a function of parallax. Note that the parallax axis is
inverted so that more distant stars are at the top. The horizontal lines show
the pulsar parallax and its uncertainties. Absorption between the vertical
lines at $2\pm 10$ km/s is viewed as being due to the ISM for our standard
model.  Figure 3: Background stars with both Na and Ca high velocity
absorption features (red squares), only one of Na or Ca high velocity
absorption features (red diamonds) or no high velocity absorption features
(green circles) superposed on the dual H$\alpha$/O[III] emission line image of
S147. The large cyan circles are the three $1^{\circ}$ diameter Hectospec
fields and the small yellow circles are foreground stars.  Figure 4:
Background stars with both Na and Ca high velocity absorption features (red
squares), only one of Na or Ca high velocity absorption features (red
diamonds) or no high velocity absorption features (green circles) superposed
on the 1.4 GHz Canadian Galactic Plane Survey map of S147. The large cyan
circles are the three $1^{\circ}$ diameter Hectospec fields and the small
yellow circles are foreground stars. Stars too close to the radio point
sources are not shown. The scale of this image is the same as in Fig. 3.

##  2 Sample and Observations

We selected stars from Gaia DR3 (Gaia Collaboration et al. 2016, Gaia
Collaboration et al. 2023) within $1.7^{\circ}$ of the center of the SNR (RA
$84.75^{\circ}$, Dec $27.83^{\circ}$, Green 2009). They were also selected to
have $G<15$ mag and $\varpi>0.2$ mas. We obtained estimates for the extinction
towards each star using the three dimensional combined19 mwdust models (Bovy
et al. 2016) which are based on the Drimmel et al. (2003), Marshall et al.
(2006) and Green et al. (2019) extinction distribution models. Fig. 2 shows
the distribution of the stars in extinction-corrected absolute magnitude and
color simply using the inverse parallax for the distance. Superposed are Solar
metallicity PARSEC (Bressan et al. 2012, Marigo et al. 2013, Pastorelli et al.
2020) isochrones with ages of $10^{7}$ to $10^{9}$ years in steps of $0.5$ dex
colored red for temperatures above $10^{4}$ K. Also marked are the candidate
unbound companion to the neutron star, HD 37424, and two stars, HD 36665 and
HD 37318, previously known to show high velocity Ca/Na absorption lines (Silk
& Wallerstein 1973, Phillips et al. 1981, Phillips & Gondhalekar 1983, Sallmen
& Welsh 2004). Dinçel et al. (2015) were concerned that these stars had
distance estimates that could put them in the foreground of the SNR, but their
Gaia parallaxes put them at or beyond the SNR distance (see Kochanek 2021).

Stars hotter than roughly $10^{4}$ K (i.e., O and B stars) do not have Ca II
or Na I absorption (the coolest B stars can). These are the optimal targets to
search for absorption lines from the interstellar medium (ISM) and the SNR. We
selected as our primary targets stars with $M_{G}<1.5$ mag and
$B_{P}-R_{P}<0.5$ mag, as shown in the red box in Fig. 2. We also required
total proper motions $\mu_{tot}<10$ mas/year. Higher observational priorities
were given to stars with parallaxes $0.4<\varpi<2.0$ mas, as these should most
strongly bracket the expected distance. Objects outside the box were included
at still lower priorities. We divided the targets into Bright and Faint
samples at $G<11.5$ mag to help obtain similar counts over the large
brightness range. The Bright sample was filled in with lower priority $G>13$
mag stars just to avoid having empty fibers.

We used the three pointings shown in Figs. 3 and 4. One was centered on the 3°
diameter SNR (“Center”) and the other two were centered contiguously to the
southeast (“SE Edge”) and northeast (“NE Edge”) of the center to capture the
edges, where we might expect lower velocities but higher projected column
densities. Each was observed in two spectrograph configurations to cover both
the Ca II and Na I wavelength ranges. We obtained a total of eleven
observations. We obtained all 4 basic observations (Bright/Faint, Ca II/Na I)
for the Center and SE Edge fields, plus one extra Bright, Na I observation of
the Center field. We only obtained the Bright Ca II/Na I observations for the
NE Edge field. For the Faint observations, a star with $G\simeq 12.5$ mag had
a signal-to-noise (S/N) ratio of $\simeq 60$ per pixel or $\simeq 300$ per
resolution element. For the shorter exposure Bright fields, the same S/N was
achieved for stars with $G\simeq 11$ mag. Table 1 summarizes the observations,
where $N_{spectra}$ is the number of stellar spectra obtained for each
pointing.

The observations were processed using the standard IDL pipeline for
Hectochelle data developed by the SAO Telescope Data Center
(TDC111https://lweb.cfa.harvard.edu/mmti/hectospec/hecto_pipe_report.pdf).
Briefly, the separate CCD images are combined with cosmic-ray elimination and
the target traces are extracted to form 1D spectra. Wavelength calibration is
supplied by ThAr exposures taken at the time of observation, and sky
subtraction is performed using contemporaneous fiber spectra of blank regions
within the 1 degree field. The wavelength stability of Hectochelle has been
measured to be $0.78\pm 0.78$ $\rm km~{}s^{-1}$over 3 years, with an rms
variation among fibers for a single pointing of just 0.17 $\rm km~{}s^{-1}$
(P. Cargile, priv. comm.). In total , we obtained spectra of 552 stars, of
which 259 had $M_{G}<1.5$ mag and $B_{P}-R_{P}<0.5$ mag (the red box region of
Fig. 2) and 293 were either less luminous or redder. Of the 259 luminous blue
stars, 172 (194) were good spectra showing clean Ca II (Na I) absorption
lines, and 101 (60) of these stars had high velocity absorption lines
(heliocentric velocities outside $2\pm 10$ km/s). The velocities of these
clean absorption lines are shown in Fig. 2 and listed in Table 2.

Table 1: Field information  Field/ | RA | Dec | Line | Time | $N_{spectra}$ | G mag | UT Date
---|---|---|---|---|---|---|---
Config | (J2000) |  | (s) |  | range |
Center/Faint | 05:39:08.4 | +27:48:06.2 | Ca II | 9900 | 201 | 11.5-15 | 2023-09-26
Center/Faint | ” | ” | Na I | 2160 | 203 | 11.5-15 | 2023-09-25
Center/Faint | ” | ” | Na I | 2160 | 201 | 11.5-15 | 2023-11-09
Center/Bright | ” | ” | Ca II | 2160 | 46 | 8.5-13 | 2023-11-04
Center/Bright | ” | ” | Na I | 900 | 46 | 8.5-13 | 2023-11-04
SE Edge/Faint | 05:42:55.9 | +27:14:48.5 | Ca II | 9360 | 203 | 11.5-15 | 2023-09-27
SE Edge/Faint | ” | ” | Na I | 3600 | 203 | 11.5-15 | 2023-10-12
SE Edge/Bright | ” | ” | Na I | 900 | 56 | 8.5-13 | 2023-11-04
SE Edge/Bright | ” | ” | Ca II | 2160 | 56 | 8.5-13 | 2023-11-04
NE Edge/Bright | 05:42:49.1 | +28:27:27.9 | Na I | 1080 | 46 | 8.5-13 | 2023-11-08
NE Edge/Bright | ” | ” | Ca II | 2160 | 46 | 8.5-13 | 2023-11-24

The spectra were individually examined to identify the velocities of Ca II and
Na I absorption features, which were flagged as being clean or not. Since
equivalent widths, particularly with the complications of overlapping, weaker
stellar and interstellar features, are unimportant for our present objectives,
we simply measured line centroids. The Ca II absorption lines were generally
stronger than the corresponding Na I absorption lines. The biggest distinction
is simply between the luminous, blue, primary targets and the fainter or
redder low priority and filler targets. Fig. 6 shows the spectra of a primary
target near the lower (reddest), right (faintest) corner of the red selection
box in Fig. 2. The star has strong interstellar absorption near zero velocity
and then strong Ca II lines near $-86$, $-69$ and $+85$ km/s. Only the $-86$
km/s line appears in the corresponding Na I spectrum. Fig. 6 shows the spectra
of a red giant star a little below the red clump in Fig 2. The low blue flux
of the red giant leads to a very noisy Ca II spectrum. The Na I spectrum has a
high signal-to-noise ratio, but in addition to the interstellar Na I
absorption, there is stellar Na I absorption plus additional stellar
atmospheric absorption lines of other species. For these stars, high velocity
interstellar absorption can only be identified after modelling the spectrum to
identify and remove these stellar features. For this paper we only use the
luminous blue stars with clean absorption features.

Figure 5: (Left) The Ca II and Na I spectra of the luminous, blue star near
the lower, right corner of the selection box in Fig. 2. This star has strong
ISM absorption near zero velocity for both species. There are three additional
strong Ca II absorption features at $-86$, $-69$, and $+85$ km/s, while only
the feature at $-86$ km/s is seen in the Na I spectrum. The strong absorption
line at $-300$ km/s is the other component of the interstellar Na I doublet at
$0$ km/s.  Figure 6: (Right) The Ca II and Na I spectra of the cool star just
below the red clump in Fig 2. There are few counts in the Ca II spectrum
because the star is a red giant. The high-signal-to-noise ratio Na I spectrum
has strong interstellar absorption near $0$ km/s, but also strong stellar
absorption lines at $+59$ km/s of multiple species including Na I, as labeled
and marked by the dashed lines and asterisks.

Figure 7: (Left) Normalized $P/P_{max}$ probability distributions for the
distance to S147 using only the Ca II (CaII, black dashed), only the Na I
(NaI, black dotted) or both (CaII$+$NaI, red solid) absorption lines for an
ISM velocity window of $2\pm 10$ km/s. The black solid curve is the distance
probability distribution implied by the pulsar parallax. The points show the
Bailer-Jones et al. (2021) photometric (upper) and photogeometric (lower)
distances to the candidate unbound binary companion HD 37424.  Figure 8:
(Right) Absorption velocities relative to the LSR as a function of offset from
the center of the SNR. The curves show the maximum expected velocities (Eqn.
4) for a radius of $R_{SNR}=1.5^{\circ}$ and velocities at the edge of
$v_{0}=100$, $150$ and $200$ km/s. The median 2 km/s LSR velocity has been
subtracted.

Figure 9: (Left) Integral distributions of the background stars with high
velocity absorption (solid red) and without high velocity absorption (dashed
black) in normalized counts from the emission line image in Fig. 3. The dotted
curve shows the distribution of the foreground stars. The Kolmogorov-Smirnov
test probability that the two background distributions are the same is
$1.5\%.$ Figure 10: (Right) Integral distributions of the background stars
with high velocity absorption (solid red) and without high velocity absorption
(dashed black) in normalized flux from the 1.4 GHz image in Fig. 4. The dotted
curve shows the distribution of the foreground stars. The Kolmogorov-Smirnov
test probability that the two background distributions are the same is 53%.

##  3 Results

Fig. 2 shows the velocities of the clean absorption components in the
luminous, blue stars as a function of parallax. There is a tight band at low
velocities due to the ISM. Because S147 lies almost directly towards the
Galactic anti-center, there is no systematic distance dependence to the ISM
velocity relative to the Sun. If we roughly exclude the higher velocity
absorption lines, the median of the ISM velocity band is $2$ km/s. For our
standard models we will assign absorption lines within $2\pm 10$ km/s to the
ISM, although we report the consequences of different choices for the velocity
width used to define the interstellar absorption There are some higher
velocity features for stars closer than the pulsar of unknown origin, but
there is a clear change in the numbers of high velocity absorption components
at the parallax of the pulsar, with velocities approaching $\sim 150$ km/s.
Kirshner & Arnold (1979) and Ren et al. (2018) found a similar velocity range
range for the emission lines from the SNR.

The goal is to determine the distance $d_{0}$ (in kpc) to the SNR, but we fit
the measured parallaxes of the stars using $\varpi_{0}=1/d_{0}$. The property
of the SNR is that it produces a change in the statistics of high velocity
absorption lines, so we model the probability of a star with distance
$d>d_{0}$ ($d<d_{0}$) having a high velocity absorption line as $p_{far}$
($p_{near}$). For a given distance, there are $n_{h}$ ($m_{h}$) $d>d_{0}$
($d<d_{0}$) stars with a high velocity absorption and $n_{l}$ ($m_{l}$)
without a high velocity feature. The probability of the observed distribution
is then

$P(D|d_{0},p_{far},p_{near})=p_{far}^{n_{h}}\left(1-p_{far}\right)^{n_{l}}p_{near}^{m_{h}}\left(1-p_{near}\right)^{m_{l}}.$
(1)

Arguably, we should include the statistical weights so that

$P(D|d_{0},p_{far},p_{near})={(n_{h}+n_{l})!(m_{h}+m_{l})!\over
n_{h}!n_{l}!m_{h}!m_{l}!}p_{far}^{n_{h}}\left(1-p_{far}\right)^{n_{l}}p_{near}^{m_{h}}\left(1-p_{near}\right)^{m_{l}},$
(2)

and we comment on the consequences of doing so below. We assume uniform priors
of $0\leq d_{0}\leq 10$ kpc, $0\leq p_{far}\leq 1$ and $0\leq p_{near}\leq 1$
and use Bayesian methods to compute the final likelihood distributions. We
include one additional elaboration to incorporate parallax uncertainties. For
a star with parallax $\varpi_{i}\pm\sigma_{i}$ and a high velocity absorption
line, we make its contribution to (for example) $n_{h}$ be

$\delta
n_{h}={1\over\sqrt{2\pi}\sigma_{i}}\int_{0}^{\varpi_{0}}d\varpi\exp\left(-{1\over
2}\left({\varpi-\varpi_{i}\over\sigma_{i}}\right)^{2}\right)$ (3)

and its contribution to $m_{h}$ would then be $1-\delta n_{h}$. This has the
particular effect of down weighting the importance of stars with large
parallax errors.

Fig. 8 shows the resulting estimates for the distance after marginalizing over
the probabilities $p_{far}$ and $p_{near}$ of having a high velocity
absorption line. The results using only the Ca II, only the Na I or both sets
of absorption features give consistent results. The Na I results are
significantly worse simply because fewer high velocity Na I features were
detected. The median distance estimates for Ca-only, Na-only and both are
$1.34$, $1.47$, and $1.37$ kpc with 90% confidence regions of $1.28$-$1.43$,
$1.18$-$1.94$ and $1.30$-$1.47$ kpc, respectively. If we increase the width of
the velocity window defining ISM absorption to $2\pm 15$ km/s or $2\pm 20$
km/s, the median distances are both $1.36$ kpc with 90% confidence ranges of
$1.29$-$1.45$ kpc and $1.28$-$1.46$ kpc, respectively, so the distance
estimate is robust to changes in the definition of low and high velocities.
Similarly, if we use Eqn. 2 instead of Eqn. 1, the median distances combining
both lines for the three velocity ranges are $1.37$, $1.35$ and $1.36$ kpc,
respectively, with 90% confidence ranges of $1.31$-$1.43$, $1.30$-$1.41$ and
$1.30$-$1.42$ kpc, so this choices matters only at the level of $\sim 0.01$
kpc. For the parallax of the pulsar ($\varpi=0.72\pm 0.12$ mas), the median
distance is 1.46 kpc with a 90% range of $1.12$-$2.10$ kpc

Fig. 8 shows the absorption line velocities relative to the LSR as a function
of each star’s offset $R$ from the center of the SNR. If the absorbing
material is in spherical, homologous expansion, then the velocities should be
confined within the curve

$v=v_{0}\left(1-{R^{2}\over R_{SNR}^{2}}\right)^{1/2}$ (4)

where $v_{0}$ is the expansion velocity at $R_{SNR}$. Such curves are shown in
Fig. 8 for $v_{0}=100$, $150$ and $200$ km/s and $R_{SNR}=1.5^{\circ}$ ($34$
pc given our median distance estimate). If we adopt $v_{0}=150$ km/s and
assume Sedov-Taylor expansion, then the age of the remnant is $94.6\pm 13.4$
thousand years with $\log(E_{51}/n_{1})=0.36\pm 0.16$, where $E_{51}$ is the
explosion energy in units of $10^{51}$ erg, $n_{0}$ is the ISM density in
units of cm-3 and we have assumed 10% errors on the radius and velocity. These
are similar to the parameters proposed from the eROSITA X-ray observations of
S146 (Michailidis et al. 2024). This age is moderately longer than the
kinematic age estimated using the pulsar distance from the SNR center and its
proper motion and significantly shorter than the pulsar spin down age
($33_{-9}^{+17}$ and 618 thousand years, respectively, Kramer et al. 2003).
The low ionization atomic Ca II and Na I states are expected only in the
colder, outermost shell of the SNR and the complex emission line morphology in
Fig. 3 suggests that the observed velocities may be subject to projection
effects where we see a higher absorption column density along lines of sight
perpendicular to shock fronts. This will tend to bias the observed absorption
velocities to be lower than for a purely spherical expansion.

Finally, we looked for correlations between the presence of high velocity
absorption lines and the SNR radio and emission line fluxes. Unfortunately,
there is no sufficiently deep, publicly available X-ray image of S147 (the
eROSITA image of Michailidis et al. 2024 is probably suitable but it is not
public). For H$\alpha$ we first tried the map from Finkbeiner (2003) which has
the limitation that its resolution is only 6’. We were serendipitiously
contacted by an astrophotographer (C. Mauche) who had just generated an
H$\alpha$ and O [III] image of S147 using an ZWO ASI29MC Pro camera and an
HUTECH/IDAS NBZ UHS dual narrow-band filter. The stars are removed by an
artificial intelligence application in StarNetV2 and the image was binned to
18″. We used the 1′ resolution 1.4 GHz map from the Canadian Galactic Plane
Survey (CGPS, Taylor et al. 2003). Figs. 3 and 4 show the dual emission line
and radio images with the three Hectospec fields and the target stars coded by
whether they are foreground, background stars with high velocity absorption or
background stars without high velocity absorption.

Figs. 10 and 10 show the normalized emission line and radio flux distributions
for each class of stars. The distributions in emission line flux for the
background stars with and without high velocity absorption are strikingly
different - stars with high velocity absorption are correlated with brighter
line emission. The distribution of the foreground stars then lies between the
two background star distributions. The Kolmogorov-Smirnov test likelihood that
the two background distributions are drawn from the same distribution is 1.5%.
The emission line image in Fig. 3 was obtained as an astrophotography project
and so is not a clean H$\alpha$ or O [III] image. However, if we use the well-
defined, but much lower resolution Finkbeiner (2003) H$\alpha$ image, we find
the same result. The distributions are similar to Fig. 10 but the Kolmogorov-
Smirnov significance of the differences is 21%, almost certainly due to the
much lower resolution. For the normalized radio flux distribution in Fig. 10,
we removed stars too close to the numerous point sources seen in Fig. 4. While
the resulting flux distributions in Fig. 10 show the same basic ordering as
the emission line flux distributions in Fig. 10, the differences are far
smaller and the Kolmogorov-Smirnov significance is only 53%.

##  4 Discussion

We believe that the idea from Cha et al. (1999) for the Vela SNR that the
appearance of high velocity Ca II and Na I in stars as a function distance can
be used to determine the distance to SNRs is now well-validated by this work.
We obtain a distance to S147 of $1.37$ kpc ($1.30$-$1.47$ at 90% confidence),
consistent with the less precise pulsar parallax and the parallax of the
candidate unbound binary companion to the SN progenitor. As seen in Fig. 3, we
obtained only one partial and two complete Hectospec pointings covering only
$\sim 20\%$ of the remnant, so it would be easy to improve the statistics
(Fig. 2 and 8), to improve the kinematics of the SNR (Fig. 8), or to better
study correlations between the absorption and other properties of the SNR
(Figs. 10 and 10).

One of our goals was to determine which absorption lines were better for this
method, and it appears to clearly be the Ca II lines. Of the 259 luminous hot
stars observed ($M_{G}<1.5$ mag, $B_{P}-R_{P}<0.4$ mag), 101 showed high
velocity Ca II absorption lines and only 60 showed high velocity Na I
absorption. The Ca II absorption velocities are also systematically higher,
with a median velocity after excluding the ISM lines of $58$ km/s relative to
the LSR, versus $25$ km/s for Na I. This is a result of needing higher total
column densities to detect Na I than Ca II, because Na I is easily
photoionized by near UV photons below 2400 Å, while Ca II requires photons
shortward of 1045 Å. Ca II and Na I originate in radiative shock waves in the
outer layers of the SNR (Cox, 1972a). Projection effects lead to lower column
densities and higher velocities normal to the shock front and the reverse
parallel to the shock front. This means that an absorption feature requiring a
higher column density will be biased towards lower average velocities. While
Ca II shows more and stronger absorption, it has the disadvantage of a bluer
wavelength (3934/3968Å versus 5890/5896Å for Na I), where extinction losses
are higher and spectrographs/detectors tend to be less efficient.

We focused on luminous, hot stars because they lack stellar Na I and Ca II
absorption (and narrow stellar absorption features generally) due to their
temperatures. This greatly simplifies identifying absorption lines from the
interstellar medium and the SNR. It is possible to apply this method to Na I
absorption in the more numerous, luminous red stars, but it is complicated by
the stellar absorption features of Na I and many other species, as seen in
Fig. 6. Essentially, the stellar spectrum has to be modelled to confidently
distinguish interstellar and stellar absorption features. This is almost
certainly feasible, but beyond our present scope. Being able to apply this
method to luminous red giants will be increasingly important for SNRs with
higher extinctions.

In old SNRs such as S147, we expect strong correlations between the Ca II and
Na I column and the H$\alpha$ and radio fluxes. The radiative shock waves
responsible for the Ca II and Na I also drive the H$\alpha$ emission (Cox,
1972b), and the brightest radio flux in old SNRs originates in radiative
shocks where the cooling, post-shock gas compresses both the ambient cosmic
ray electrons and the magnetic field (van der Laan, 1962; Raymond et al.,
2020a; Tutone et al., 2021). We correlated the background stars with and
without high velocity absorption lines with emission line and radio images of
S147. We found a significant correlation of high velocity absorption with the
flux in an H$\alpha$/O[III] astrophotography image of S147 and, more weakly,
with the far lower resolution H$\alpha$ image from Finkbeiner (2003). This
suggests that using emission line maps to prioritize targets is likely more
important than our center versus edge experiment. We found no significant
correlation between the background stars with high velocity absorption lines
and the Canadian Galactic Plane Survey (Taylor et al. 2003) 1.4 GHz radio map
of S147. The insignificant correlation with the radio emission may be driven
by a combination of the modest (1′) resolution and limited dynamic range of
the CGPS radio map for the emission from the SNR. Unfortunately, no
sufficiently deep, publicly available X-ray image of S147 was available to
check for correlations.

Based on these results, we initiated a program to obtain data for up to
sixteen additional Northern and Southern SNRs using MMT/Hectospec and
Magellan/M2FS in the first trimester of 2024. As emphasized in the
introduction, modern multi-fiber, echelle spectrographs allow these
observations to be carried out very efficiently. The selection of target stars
for this trimester is not making use of the correlation between the existence
of high velocity absorption and emission line flux we find here, but will
certainly serve to confirm the correlation for use in future observations. We
are also investigating modeling the spectra of the red giants, which are
currently just used to fill otherwise unused fibers, to also use them in the
search for absorption features from the SNRs.

We thank Christopher Mauche for sharing his images of S147 with us. We also
had helpful discussions regarding sodium absorption with Jesse Han. CSK is
supported by NSF grants AST-1908570 and AST-2307385. Observations reported
here were obtained at the MMT Observatory, a joint facility of the Smithsonian
Institution and the University of Arizona. Data was reduced by the SAO TDC; we
thank Sean Moran and Jaehyon Rhee for their help in that work.

## References

  * Badenes et al. (2009) Badenes, C., Harris, J., Zaritsky, D., & Prieto, J. L. 2009, ApJ, 700, 727, doi: 10.1088/0004-637X/700/1/727
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Blundell & Bowler (2004) Blundell, K. M., & Bowler, M. G. 2004, ApJ, 616, L159, doi: 10.1086/426542
  * Bovy et al. (2016) Bovy, J., Rix, H.-W., Green, G. M., Schlafly, E. F., & Finkbeiner, D. P. 2016, ApJ, 818, 130, doi: 10.3847/0004-637X/818/2/130
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Cha et al. (1999) Cha, A. N., Sembach, K. R., & Danks, A. C. 1999, ApJ, 515, L25, doi: 10.1086/311968
  * Chatterjee et al. (2009) Chatterjee, S., Brisken, W. F., Vlemmings, W. H. T., et al. 2009, ApJ, 698, 250, doi: 10.1088/0004-637X/698/1/250
  * Corbet et al. (2011) Corbet, R. H. D., Cheung, C. C., Kerr, M., et al. 2011, The Astronomer’s Telegram, 3221, 1
  * Cox (1972a) Cox, D. P. 1972a, ApJ, 178, 143, doi: 10.1086/151774
  * Cox (1972b) —. 1972b, ApJ, 178, 159, doi: 10.1086/151775
  * Dinçel et al. (2015) Dinçel, B., Neuhäuser, R., Yerli, S. K., et al. 2015, MNRAS, 448, 3196, doi: 10.1093/mnras/stv124
  * Dodson et al. (2003) Dodson, R., Legge, D., Reynolds, J. E., & McCulloch, P. M. 2003, ApJ, 596, 1137, doi: 10.1086/378089
  * Drimmel et al. (2003) Drimmel, R., Cabrera-Lavers, A., & López-Corredoira, M. 2003, A&A, 409, 205, doi: 10.1051/0004-6361:20031070
  * Fabricant et al. (2005) Fabricant, D., Fata, R., Roll, J., et al. 2005, PASP, 117, 1411, doi: 10.1086/497385
  * Fesen et al. (2018) Fesen, R. A., Neustadt, J. M. M., Black, C. S., & Milisavljevic, D. 2018, MNRAS, 475, 3996, doi: 10.1093/mnras/sty072
  * Fesen et al. (2021) Fesen, R. A., Weil, K. E., Cisneros, I., Blair, W. P., & Raymond, J. C. 2021, MNRAS, 507, 244, doi: 10.1093/mnras/stab2066
  * Finkbeiner (2003) Finkbeiner, D. P. 2003, ApJS, 146, 407, doi: 10.1086/374411
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Green (2009) Green, D. A. 2009, Bulletin of the Astronomical Society of India, 37, 45, doi: 10.48550/arXiv.0905.3699
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93, doi: 10.3847/1538-4357/ab5362
  * Hinton et al. (2009) Hinton, J. A., Skilton, J. L., Funk, S., et al. 2009, ApJ, 690, L101, doi: 10.1088/0004-637X/690/2/L101
  * Jennings et al. (2014) Jennings, Z. G., Williams, B. F., Murphy, J. W., et al. 2014, ApJ, 795, 170, doi: 10.1088/0004-637X/795/2/170
  * Kameswara Rao et al. (2020) Kameswara Rao, N., Lambert, D. L., Reddy, A. B. S., et al. 2020, MNRAS, 493, 497, doi: 10.1093/mnras/staa231
  * Kerzendorf et al. (2009) Kerzendorf, W. E., Schmidt, B. P., Asplund, M., et al. 2009, ApJ, 701, 1665, doi: 10.1088/0004-637X/701/2/1665
  * Kirshner & Arnold (1979) Kirshner, R. P., & Arnold, C. N. 1979, ApJ, 229, 147, doi: 10.1086/156939
  * Kochanek (2021) Kochanek, C. S. 2021, MNRAS, 507, 5832, doi: 10.1093/mnras/stab2483
  * Kochanek (2022) —. 2022, MNRAS, 511, 3428, doi: 10.1093/mnras/stac098
  * Kochanek et al. (2019) Kochanek, C. S., Auchettl, K., & Belczynski, K. 2019, MNRAS, 485, 5394, doi: 10.1093/mnras/stz717
  * Kramer et al. (2003) Kramer, M., Lyne, A. G., Hobbs, G., et al. 2003, ApJ, 593, L31, doi: 10.1086/378082
  * Lin et al. (2023) Lin, R., van Kerkwijk, M. H., Kirsten, F., Pen, U.-L., & Deller, A. T. 2023, ApJ, 952, 161, doi: 10.3847/1538-4357/acdc98
  * Margon (1984) Margon, B. 1984, ARA&A, 22, 507, doi: 10.1146/annurev.aa.22.090184.002451
  * Marigo et al. (2013) Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488, doi: 10.1093/mnras/stt1034
  * Marshall et al. (2006) Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635, doi: 10.1051/0004-6361:20053842
  * Mateo et al. (2012) Mateo, M., Bailey, J. I., Crane, J., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84464Y, doi: 10.1117/12.926448
  * Michailidis et al. (2024) Michailidis, M., Pühlhofer, G., Becker, W., et al. 2024, arXiv e-prints, arXiv:2401.17312, doi: 10.48550/arXiv.2401.17312
  * Ng et al. (2007) Ng, C. Y., Romani, R. W., Brisken, W. F., Chatterjee, S., & Kramer, M. 2007, ApJ, 654, 487, doi: 10.1086/510576
  * Pasquini et al. (2000) Pasquini, L., Avila, G., Allaert, E., et al. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 129–140, doi: 10.1117/12.395491
  * Pastorelli et al. (2020) Pastorelli, G., Marigo, P., Girardi, L., et al. 2020, MNRAS, 498, 3283, doi: 10.1093/mnras/staa2565
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49
  * Phillips & Gondhalekar (1983) Phillips, A. P., & Gondhalekar, P. M. 1983, MNRAS, 202, 483, doi: 10.1093/mnras/202.2.483
  * Phillips et al. (1981) Phillips, A. P., Gondhalekar, P. M., & Blades, J. C. 1981, MNRAS, 195, 485, doi: 10.1093/mnras/195.3.485
  * Raymond et al. (2020a) Raymond, J. C., Chilingarian, I. V., Blair, W. P., et al. 2020a, ApJ, 894, 108, doi: 10.3847/1538-4357/ab886d
  * Raymond et al. (2020b) Raymond, J. C., Slavin, J. D., Blair, W. P., et al. 2020b, ApJ, 903, 2, doi: 10.3847/1538-4357/abb821
  * Ren et al. (2018) Ren, J.-J., Liu, X.-W., Chen, B.-Q., et al. 2018, Research in Astronomy and Astrophysics, 18, 111, doi: 10.1088/1674-4527/18/9/111
  * Ritchey et al. (2024) Ritchey, A. M., Federman, S. R., & Lambert, D. L. 2024, MNRAS, 528, 4490, doi: 10.1093/mnras/stae180
  * Sallmen & Welsh (2004) Sallmen, S., & Welsh, B. Y. 2004, A&A, 426, 555, doi: 10.1051/0004-6361:20040093
  * Sheinis (2016) Sheinis, A. I. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 99081C, doi: 10.1117/12.2234334
  * Silk & Wallerstein (1973) Silk, J., & Wallerstein, G. 1973, ApJ, 181, 799, doi: 10.1086/152090
  * Szentgyorgyi et al. (2011) Szentgyorgyi, A., Furesz, G., Cheimets, P., et al. 2011, PASP, 123, 1188, doi: 10.1086/662209
  * Taylor et al. (2003) Taylor, A. R., Gibson, S. J., Peracaula, M., et al. 2003, AJ, 125, 3145, doi: 10.1086/375301
  * Temim et al. (2022) Temim, T., Slane, P., Raymond, J. C., et al. 2022, ApJ, 932, 26, doi: 10.3847/1538-4357/ac6bf4
  * Tutone et al. (2021) Tutone, A., Ballet, J., Acero, F., D’Aì, A., & Cusumano, G. 2021, A&A, 656, A139, doi: 10.1051/0004-6361/202141978
  * van der Laan (1962) van der Laan, H. 1962, MNRAS, 124, 125, doi: 10.1093/mnras/124.2.125

Table 2: Absorption Velocities of the Hot, Blue Stars

Gaia DR3 ID | Ca II | Na I | Gaia DR3 ID | Ca II | Na I
---|---|---|---|---|---
| (km/s) | (km/s) |  | (km/s) | (km/s)
3441702193599182336 |  $-118$,$4$,$38$ | $2$ | 3441703598051819776 | $3$ | $2$
3441706518633166080 |  |  $-21$,$0$ | 3441723496636807040 |  |
3441698860704616064 |  |  | 3442475974907331200 | $3$ | $0$
3442474561861744000 |  |  $-13$,$0$ | 3442449552266565888 | $3$ | $5$
3442476112346300800 | $3$ | $2$ | 3442452197968495616 |  $-105$,$5$ | $3$
3442455011170611712 | $5$ | $5$ | 3442479337865371392 |  $-32$,$1$,$87$ |  $-34$,$0$
3442469751498370560 |  $-127$,$2$ | $3$ | 3442472779451400448 | $3$ | $0$
3442491947889351296 | $3$ | $4$ | 3442491814744678528 |  | $5$
3442546824687824256 | $0$ | $3$ | 3442506348915651328 |  $5$,$23$ | $2$
3442599669964938368 |  $-58$,$0$ |  $-14$,$-2$ | 3442503527120074624 |  |  $-24$,$1$

Note. — The entries are the velocities of the Ca II and Na I absorption lines
separated by commas. The complete table is available in the electronic version
of the paper.


# Necessary and sufficient stability conditions for neutral-type delay
systems: Polynomial approximations

Gerson Portilla111The work of the first and third authors was supported by
project CONACYT A1-S-24796, Mexico<EMAIL_ADDRESS>Mathieu
Bajodek<EMAIL_ADDRESS>Sabine Mondié222The work of the first and
third authors was supported by project CONACYT A1-S-24796, Mexico.
<EMAIL_ADDRESS>Department of Automatic Control, CINVESTAV-IPN,
07360 Mexico D.F., Mexico  Laboratoire d’Automatique, de Génie des Procédés et
de Génie Pharmaceutique (LAGEPP), CPE Lyon, UCBL, Campus Lyontech-La Doua,
69100 Villeurbanne, France

###### Abstract

A new necessary and sufficient stability test in a tractable number of
operations for linear neutral-type delay systems is introduced. It is
developed in the Lyapunov-Krasovskii framework via functionals with prescribed
derivative. The necessary conditions, which stem from substituting any
polynomial approximation of the functional argument, reduce to a quadratic
form of monomials whose matrix is independent of the coefficients of the
approximation under consideration. In the particular case of Chebyshev
polynomials, the functional approximation error is quantified, leading to an
estimate of the order of approximation such that the positive semi-
definiteness of the functional is verified. Some examples illustrate the
obtained results.

###### keywords:

Neutral type delay systems, stability criterion, delay Lyapunov matrix.

††journal: Journal of the Franklin Institute

##  1 Introduction

In the Lyapunov-Krasovskii framework, stability theorems for neutral-type
systems were introduced by Krasovskii in [1]. In the monograph by Kharitonov
[2], functionals with prescribed derivative for linear neutral time-delay
systems are constructed. They rely on the so-called delay Lyapunov matrix
that, in analogy with delay-free linear systems, is obtained as the solution
of a Lyapunov equation, which is now a set of three properties called
symmetry, dynamic and algebraic. These functionals with prescribed derivative
have been exploited to present sufficient [3], necessary and sufficient
stability [4, 5] results for neutral type delay systems. The above-mentioned
results allowed the presentation of stability and instability theorems on a
special set of bounded functional arguments as in [3], which turned out to be
crucial for achieving tractable sufficiency results.

Necessary conditions follow from the substitution of the functional argument.
In particular, substituting the functional argument by an approximation in
terms of the system’s fundamental matrix reveals an elegant necessary
condition in terms of the delay Lyapunov matrix [4].

For assessing the sufficiency of these conditions, one must find the order of
approximation of an arbitrary initial function that satisfies the positivity
condition of the functional for any initial approximated function. The
fundamental matrix approximation is poor, thus it leads to very large
approximation orders guaranteeing sufficiency. To remedy this issue, authors
have proposed, to the cost of losing the delay Lyapunov matrix formulated
criteria, better approximations of the functional argument, for example, the
piece-wise linear approximations [6] and Legendre polynomials in [7] proposed
for the retarded type case. The orthogonality, speed of convergence, and
accuracy of Legendre polynomial approximations allowed achieving similar or
tighter approximation orders piece-wise linear approximations [6]. In these
contributions, after calculating the delay Lyapunov matrix, its derivative,
and substitution of the approximated argument, the computation of integrals
leads to an approximation of the functional, which, combined with sufficiency
stability theorems, gives tractable conditions in the form of linear matrix
inequalities [8, 9].

Here, considering that it is well-known that polynomial approximations perform
well and that any functional argument allows presenting necessary stability
conditions for linear neutral time-delay systems, we propose an approach
encompassing all classes of polynomial approximations. The result relies on
the fact that any polynomial approximation can be reorganized as the product
of a monomial vector multiplied by a vector of appropriate coefficients. The
positive semi-definiteness of the resulting quadratic form delivers a
necessary stability test for general polynomial approximation-based stability
results. It is characterized by integrals of the delay Lyapunov matrix
multiplied by monomials, which can be computed by a recursive method following
the one introduced in [6].

For sufficiency, as the argument approximation order influence on the
functional error must be quantified, one must restrict the analysis to a
particular approximation method. Among polynomial approximations [10], Taylor
expansion, Bessel series, Lagrange interpolations, Padé approximants, and
orthogonal polynomial approximations, such as Chebyshev and Legendre, are
available. We choose Chebyshev polynomials because, while they share the
orthogonal property with Legendre polynomials, they outperform them with the
super-geometric convergence of their approximation coefficients, ensuring a
fast error convergence rate and simplification of the computations [11]. Some
examples allow evidencing the reduced approximation orders, which are compared
with previous results in the literature for neutral-type systems.

This paper is organized as follows. In Section 2, some preliminaries on linear
neutral time-delay systems, functionals with prescribed derivative, and an
instrumental instability result are introduced. Also, a monomial
representation for any polynomial approximation and convergence properties of
Chebyshev polynomials are presented. In Section 3, a new necessary stability
test resulting from general polynomial approximations and a stability
criterion in a tractable number of operations is obtained for the case of
Chebyshev polynomials. At the end, the technical issue of computing the
integrals in terms of the delay Lyapunov matrix is tackled via a new recursive
method described in Section 4 and some examples in Section 5 validate the
results.

Notation:  We consider the spaces of $\mathbb{R}^{n}$-valued piecewise
continuously differentiable and smooth functions on $[-h,0]$, which are
denoted by $PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ and
$C_{\infty}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)},$ respectively. They are
equipped with the uniform norm

$\|\varphi\|_{h}=\sup_{\theta\in[-h,0]}\|\varphi(\theta)\|,$

where $\|\cdot\|$ stands for the Euclidean norm for vectors and the spectral
norm for matrices. $\Re(s)$ denotes the real part of a complex value $s$;
$\lambda_{\textup{min}}(W)$ is the smallest eigenvalue of a square matrix $W$;
notation $k=\overline{n_{1},n_{2}}$, where $n_{1},~{}n_{2}\in\mathbb{Z}$,
$n_{1}<n_{2}$, means that $k$ is an integer between $n_{1}$ and $n_{2}$;
$I_{n}$ stands for the $n\times n$ identity matrix; $\mathrm{He}(M)$ means
$M+M^{\top}$; $\lceil\cdot\rceil$ denotes the ceiling function;
$\\{A_{ij}\\}_{i,j=1}^{r}$ denotes a square block matrix, where
$A_{ij}\in\mathbb{R}^{n\times n},~{}i,j=\overline{1,r}$, is the block in the
$i$-th row and the $j$-th column; $\textup{vec}(X)$ means the column
vectorization of a matrix $X$; Following [2], $A\otimes B$ stands for the
Kronecker product, namely,

$A\otimes B=\begin{pmatrix}b_{11}A&b_{21}A&\cdots&b_{n1}A\\\
b_{12}A&b_{22}A&\cdots&b_{n2}A\\\ \vdots&\vdots&\ddots&\vdots\\\
b_{1n}A&b_{2n}A&\cdots&b_{nn}A\\\ \end{pmatrix},$

where $B=\\{b_{ij}\\}_{i,j=1}^{n},~{}b_{ij}\in\mathbb{R}$. Notice that, with
this definition, $\textup{vec}(AXB)=(A\otimes B)\textup{vec}(X)$. For a
symmetric matrix $\Lambda,$ the notation $\Lambda>0$ $(\Lambda\geq 0)$ means
that $\Lambda$ is a positive definite (positive semidefinite) matrix.
$\mathcal{W}(z)$ denotes the Lambert function given by $\mathcal{W}:z\to
y,~{}z\in\mathbb{R}_{+},~{}y\in\mathbb{R}_{+}$, which is uniquely defined by
the relation $y\textup{e}^{y}=z$.

##  2 Preliminaries

###  2.1 Neutral-type delay systems and Lyapunov-Krasovskii functional

Consider a neutral-type linear time-delay system of the form

$\frac{\mathrm{d}}{\mathrm{d}t}[x(t)-Dx(t-h)]=A_{0}x(t)+A_{1}x(t-h),$ (1)

where $h>0$, $A_{0}$, $A_{1}$, and $D$ are given real $n\times n$ matrices.
For $\varphi\in PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$, the solution
$x(t)=x(t,\varphi)$ is a piecewise continuous function that satisfies system
(1) almost everywhere for $t\geqslant 0$, and the difference $x(t)-Dx(t-h)$ is
continuous for $t\geqslant 0,$ except for possibly a countable number of
points. The restriction of the solution $x(t,\varphi)$ to the interval
$[t-h,t],~{}t\geqslant 0$, is denoted by

$x_{t}(\varphi):\theta\mapsto x(t+\theta,\varphi),~{}\theta\in[-h,0].$

######  Definition 1.

System (1) is exponentially stable if there exist $\gamma>0$ and $\sigma>0$
such that for any initial function $\varphi\in
PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$,

$\|x(t,\varphi)\|\leqslant\gamma e^{-\sigma t}\|\varphi\|_{h},~{}t\geqslant
0.$

We also introduce the following assumption on matrix $D$, which is a necessary
condition for the exponential stability of system (1).

######  Assumption 1.

Matrix $D$ is a Schur stable matrix, it satisfies $\|D\|<1$.

The delay Lyapunov matrix definition and a stability result on a particular
set of initial functions in the Lyapunov-Krasovskii framework are recalled.

######  Definition 2.

[2] Let $W\in\mathbb{R}^{n\times n}$ be a positive definite matrix. The _delay
Lyapunov matrix_ $U\\!:[-h,h]\rightarrow\mathbb{R}^{n\times n}$ is a
continuous matrix function, which satisfies the following properties.

  1. 1.

Dynamic property

$U^{\prime}(\theta)-U^{\prime}(\theta-h)D=U(\theta)A_{0}+U(\theta-h)A_{1},~{}\theta\in(0,h).$
(2)

  2. 2.

Symmetry property

$U^{\top}(\theta)=U(-\theta),~{}\theta\in[-h,h].$ (3)

  3. 3.

Algebraic property

$P-D^{\top}PD=-W,$ (4)

where $P=\lim\limits_{\theta\rightarrow
0}\left(U^{\prime}(+\theta)-U^{\prime}(-\theta)\right)$.

Observe that, for negative values of the argument, the delay Lyapunov matrix
satisfies

$U^{\prime}(\theta)-D^{\top}U^{\prime}(\theta+h)=-A_{0}^{\top}U(\theta)-A_{1}^{\top}U(\theta+h),~{}\theta\in(-h,0).$
(5)

The following result gives conditions for the existence and uniqueness of the
delay Lyapunov matrix. Notice that it does not rely on the stability or
instability of system (1).

######  Theorem 1.

[2] System (1) admits a unique Lyapunov matrix if and only if the system
satisfies the Lyapunov condition, i.e., if there exists $\varepsilon>0$ such
that any two points $s_{1}$ and $s_{2}$ of the spectrum of system (1),

$\Lambda=\Big{\\{}s~{}\Big{|}~{}\textup{det}\Big{(}sI-
A_{0}-se^{-sh}D-e^{-sh}A_{1}\Big{)}=0\Big{\\}},$

satisfy $|s_{1}+s_{2}|>\varepsilon$.

In the sequel, we use the following Lyapunov matrix-based functional [2,
Chapter 6, p. 242]

$\displaystyle
v_{0}(\varphi)=(\varphi(0)-D\varphi(-h))^{\top}U(0)(\varphi(0)-D\varphi(-h))+\sum\limits_{j=1}^{6}I_{j},$
(6)

where

$\displaystyle I_{1}$
$\displaystyle=2(\varphi(0)-D\varphi(-h))^{\top}\int_{-h}^{0}U^{\top}(h+\theta)A_{1}\varphi(\theta)\mathrm{d}\theta$
$\displaystyle I_{2}$
$\displaystyle=-2(\varphi(0)-D\varphi(-h))^{\top}\int_{-h}^{0}U^{\prime\top}(h+\theta)D\varphi(\theta)\mathrm{d}\theta,$
$\displaystyle I_{3}$
$\displaystyle=\int_{-h}^{0}\int_{-h}^{0}\varphi^{\top}(\theta_{1})A_{1}^{\top}U(\theta_{1}-\theta_{2})A_{1}\varphi(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1},$
$\displaystyle I_{4}$
$\displaystyle=2\int_{-h}^{0}\int_{-h}^{0}\varphi^{\top}(\theta_{1})A_{1}^{\top}U^{\prime}(\theta_{1}-\theta_{2})D\varphi(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1},$
$\displaystyle I_{5}$
$\displaystyle=-\int_{-h}^{0}\int_{-h}^{\theta_{1}}\varphi^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\varphi(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}$
$\displaystyle-\int_{-h}^{0}\int_{\theta_{1}}^{0}\varphi^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\varphi(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1},$
$\displaystyle I_{6}$
$\displaystyle=-\int_{-h}^{0}\varphi^{\top}(\theta)D^{\top}PD\varphi(\theta)\mathrm{d}\theta.$

and where $U(\theta)$ and $P$ are defined by Definition 2 with matrix $W>0$.
Its time derivative along the solution of system (1) has a prescribed
quadratic negative value of the form

$\frac{\mathrm{d}v_{0}(x_{t})}{\mathrm{d}t}=-x^{\top}(t)Wx(t),\quad t\geq 0.$
(7)

On the one hand, it was shown in [3] that if system (1) is exponentially
stable, then functional $v_{0}(\varphi)$ admits a quadratic lower bound
considering a particular set of functions given by

$\mathcal{S}=\Bigl{\\{}\varphi\in
C_{\infty}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}\Bigl{\arrowvert}\|\varphi(0)\|=1,\|\varphi^{(k)}\|_{h}\leq
r^{k},\;\forall k\in\mathbb{N}\Bigr{\\}},$ (8)

where $r=\frac{\|A_{0}\|+\|A_{1}\|}{1-\|D\|}$. On the other hand, this
functional allows proving a crucial instability result reminded below.

######  Lemma 1.

[3] If system (1) is unstable, then there exists a function
$\varphi\in\mathcal{S}$ such that

$v_{0}(\varphi)<-a_{0},\quad a_{0}=\dfrac{\lambda_{\textup{min}}(W)}{4r}.$

######  Remark 1.

We can take $\Re(s)\leqslant|s|\leqslant r$ when the system is unstable, i.e.,
for any eigenvalue such that $\Re(s)>0.$

###  2.2 Preliminaries on polynomial approximation

Next, we present some results on polynomial approximation. In particular,
their monomial representation and the convergence of Chebyshev polynomial
projections.

####  2.2.1 General polynomial approximations

This section introduces a monomial representation of polynomial approximations
for the argument $\varphi\in PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ of
functional (6). This representation reorganizes and expresses any polynomial
projection as a decoupling of monomials and coefficients. Replacing this
decoupling into functional (6) leads to a quadratic form, whose matrix is
determined by monomials and the delay Lyapunov matrix (see Section 3.1).

Let us consider a function $\varphi\in
PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ and a $N$-order polynomial
approximation of this function given by

$\varphi_{N}(\theta)=P_{N}^{\top}(\theta)\mathcal{Q},\quad\theta\in[-h,0],\quad
N\in\mathbb{N},$ (9)
$P_{N}^{\top}(\theta)=\begin{bmatrix}p_{0}(\theta)I_{n}&p_{1}(\theta)I_{n}&\cdots&p_{N-1}(\theta)I_{n}\end{bmatrix},$

where $\\{p_{k}\\}_{k\in\overline{0,N-1}}$ is a set of polynomials of order
lower than $N$ and, vector $\mathcal{Q}\in\mathbb{R}^{nN}$ collects the
coefficients associated to the chosen polynomial approximation. A monomial
representation of the previous polynomial approximation can be expanded,
resulting in the following expression:

$\varphi_{N}(\theta)=\Theta_{N}^{\top}(\theta)\Phi_{N},$ (10)
$\Theta_{N}^{\top}(\theta)=\begin{bmatrix}I_{n}&\theta
I_{n}&\cdots&\theta^{N-1}I_{n}\end{bmatrix},$

where $\Theta_{N}(\theta)$ is a matrix of monomials, and
$\Phi_{N}\in\mathbb{R}^{nN}$ is a vector with suitable coefficients after
expansion. Then, it follows from (10) that any function $\varphi\in
PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ can be written as

$\begin{split}\varphi(\theta)&=\varphi_{N}(\theta)+\tilde{\varphi}_{N}(\theta),\\\
\varphi_{N}(\theta)&=\Theta_{N}^{\top}(\theta)\Phi_{N},\quad\theta\in[-h,0],\end{split}$
(11)

where the function $\varphi_{N}(\theta)$, $\theta\in[-h,0]$, denotes the
polynomial approximation of $\varphi\in
PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ given by any polynomial
approximation, and $\tilde{\varphi}_{N}(\theta)$, $\theta\in[-h,0]$, stands
for the corresponding approximation error.

######  Remark 2.

For instance, this general framework gathers Lagrange interpolation,
orthogonal projections, or Taylor series.

####  2.2.2 Orthogonal polynomial projections

In the approximation theory, the orthogonal polynomial bases are characterized
by convergence properties and accurate approximations with few terms. In this
section, we assume that the sequence $\\{p_{k}\\}_{k\in\mathbb{N}}$ in (9) is
an orthogonal polynomial basis in the Banach space
$PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ equipped with the norm
$\varphi\mapsto\int\limits_{-h}^{0}\|\varphi(\theta)\|^{2}w(\theta)\textup{d}\theta$,
where $w$ is a weight function.

We prove next a useful property of orthogonal polynomial projection on the
sequence $\\{p_{k}\\}_{k\in\mathbb{N}}$.

######  Property 1.

For any function $\varphi\in PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ and
any order $N\in\mathbb{N}$, its orthogonal polynomial projection is given by

$\varphi_{N}:\left\\{\begin{aligned} [&-h,0]\to\mathbb{R}^{n},\\\
\theta&\mapsto\varphi_{N}(\theta):=P_{N}^{\top}(\theta)D_{N}^{-1}\int_{-h}^{0}P_{N}(\tau)\varphi(\tau)w(\tau)\mathrm{d}\tau,\end{aligned}\right.$
(12)

where $D_{N}=\int_{-h}^{0}P_{N}(\tau)P_{N}^{\top}(\tau)w(\tau)\mathrm{d}\tau$.
This approximation can also be written as

$\varphi_{N}(\theta)=\Theta_{N}^{\top}(\theta)S_{N}^{-1}\int_{-h}^{0}\Theta_{N}(\tau)\varphi(\tau)w(\tau)\mathrm{d}\tau,$
(13)

where
$S_{N}=\int_{-h}^{0}\Theta_{N}(\tau)\Theta_{N}^{\top}(\tau)w(\tau)\mathrm{d}\tau.$
Here, $D_{N}$ and $S_{N}$ are symmetric Gramm-Schmidt orthonormalization
matrices.

###### Proof.

The expression (12) follows from the Galerkin-like methods introduced in [12]
and stands as a definition. Thus, we proceed to prove (13). First, we express
the monomials in terms of their orthogonal polynomials approximation at order
$N$ for any $k\in\\{0,\dots,N-1\\}$ as follows

$\theta^{k}I_{n}=P_{N}^{\top}(\theta)D_{N}^{-1}\int_{-h}^{0}P_{N}(\tau)\tau^{k}w(\tau)\mathrm{d}\tau.$

Since a polynomial approximation can reconstruct any polynomial of order $N$,
this is indeed an equality. By concatenation, we get

$\Theta_{N}^{\top}(\theta)=P_{N}^{\top}(\theta)D_{N}^{-1}\underbrace{\int_{-h}^{0}P_{N}(\tau)\Theta_{N}^{\top}(\tau)w(\tau)\mathrm{d}\tau}_{=T_{N}},$
(14)

where $T_{N}\in\mathbb{R}^{nN\times nN}$ is a nonsingular triangular matrix
that links the bases $\Theta_{N}$ and $P_{N}$. Note also that
$\left(D_{N}^{-1}\right)^{\top}=D_{N}^{-1}$ and that

$S_{N}=\int_{-h}^{0}\Theta_{N}(\tau)\Theta_{N}^{\top}(\tau)w(\tau)\mathrm{d}\tau=\underbrace{\int_{-h}^{0}\Theta_{N}(\tau)P_{N}^{\top}(\tau)w(\tau)\mathrm{d}\tau}_{T_{N}^{\top}}D_{N}^{-1}T_{N}.$

Gathering the two previous expressions yields

$\begin{split}\Theta_{N}^{\top}(\theta)S_{N}^{-1}\int_{-h}^{0}\Theta_{N}(\tau)\varphi(\tau)w(\tau)\mathrm{d}\tau&=\left(P_{N}^{\top}(\theta)D_{N}^{-1}T_{N}\right)\left(T_{N}^{-1}D_{N}(T_{N}^{\top})^{-1}\right)\\\
&\phantom{=}\times\left(T_{N}^{\top}D_{N}^{-1}\int_{-h}^{0}P_{N}(\tau)\varphi(\tau)w(\tau)\mathrm{d}\tau\right)\\\
&=P_{N}^{\top}(\theta)D_{N}^{-1}\int_{-h}^{0}P_{N}(\tau)\varphi(\tau)w(\tau)\mathrm{d}\tau.\end{split}$

Therefore, we proved that the approximation forms provided in (12) and (13)
are equivalent. ∎

Note that
$\mathcal{Q}=D_{N}^{-1}\int_{-h}^{0}P_{N}(\tau)\varphi(\tau)w(\tau)\mathrm{d}\tau$
and
$\Phi_{N}=S_{N}^{-1}\int_{-h}^{0}\Theta_{N}(\tau)\varphi(\tau)w(\tau)\mathrm{d}\tau$
for the polynomial approximations (9) and (10), respectively.

######  Remark 3.

The first expression (12) is a standard formulation that comes from orthogonal
polynomial approximation theory [11]. Due to orthogonality of the sequence
$\\{p_{k}\\}_{k\in\mathbb{N}}$, the Gramm-Schmidt orthonormalization matrix
$D_{N}$ is diagonal which simplifies the maneuverability of the approximation
structure. The second expression (13) is an the equivalent monomial
representation and is associated to the Gramm-Schmidt orthonormalization
matrix $S_{N}$. Such a structure encompasses Legendre or Chebyshev polynomial
projections.

####  2.2.3 Chebyshev orthogonal polynomial projection

Among orthogonal polynomial bases, Chebyshev and Legendre’s polynomials are
widely investigated and applied to approximation theory due to the super-
geometric convergence of their coefficients [13]. However, the Chebyshev
coefficients decay a $\sqrt{N}$ factor faster than Legendre ones [14, 15].
Therefore, we choose the Chebyshev polynomials as a suitable polynomial
approximation basis to tackle the estimation of the functional approximation
error. In this context, $\\{p_{k}\\}_{k\in\mathbb{N}}$ are the Chebyshev
polynomial, which are given by
$p_{0}(\theta)=1,~{}p_{1}(\theta)=\frac{2\theta}{h}+1,~{}p_{2}=\frac{8\theta}{h}+\frac{8\theta^{2}}{h^{2}}+1$
for low orders, and the weight is equal to
$w(\theta)=\frac{1}{\sqrt{1-(-\frac{2\theta+h}{h})^{2}}}$ for
$\theta\in[-h,0]$.

Next, a lemma on the Chebyshev convergence rate on the set of functions
$\mathcal{S}$ is stated.

######  Lemma 2.

Consider $\varphi\in\mathcal{S}$. The Chebyshev approximation error
$\tilde{\varphi}_{N}$ satisfies the following inequality

$\|\tilde{\varphi}_{N}\|_{h}\leqslant\frac{4\left(\frac{hr}{2}\right)^{N}}{N!}.$
(15)

###### Proof.

The proof is based on Theorem 4.2 in [16] with a slight modification. First of
all, we start shifting the Chebyshev polynomials and coefficients on the
interval $[-h,0]$, as follows

$a_{k}=\frac{4}{\pi
h}\int\limits_{-h}^{0}\frac{\varphi(\theta)(-1)^{k}p_{k}(-\frac{2\theta+h}{h})}{\sqrt{1-(-\frac{2\theta+h}{h})^{2}}}\mathrm{d}\theta,\quad\theta\in[-h,0].$

The change of variable $x=-\frac{2\theta+h}{h}$ yields

$a_{k}=\frac{2(-1)^{k}}{\pi}\int\limits_{-1}^{1}\frac{\varphi(-\frac{h}{2}(1+x))p_{k}(x)}{\sqrt{1-x^{2}}}\mathrm{d}x,\quad
x\in[-1,1].$

From here, one proceeds as in Theorem 4.2 in [16]. Thus, the coefficients
admit the upper bound

$|a_{k}|\leqslant\frac{2\left(\frac{hr}{2}\right)^{m+1}}{k(k-1)\cdots(k-m)},\quad
m\geqslant 1,\quad k\geqslant m+1.$

Then, we obtain the following estimate for the approximation error

$\displaystyle\lVert\tilde{\varphi}_{N}\rVert_{h}$
$\displaystyle\leqslant\sum_{k=N}^{\infty}\frac{2\left(\frac{hr}{2}\right)^{m+1}}{k(k-1)\cdots(k-m)},$
$\displaystyle=\frac{2\left(\frac{hr}{2}\right)^{m+1}}{m}\sum_{k=N}^{\infty}\left(\tfrac{1}{(k-1)\cdots(k-m)}-\tfrac{1}{k(k-1)\cdots(k-m+1)}\right),$
$\displaystyle\leqslant\frac{2\left(\frac{hr}{2}\right)^{m+1}}{m(N-1)\cdots(N-m)},\quad
m\geqslant 1,\quad N\geqslant m+1.$

Finally, by taking $m=N-1$, we conclude that for all $N\geqslant 1$,
$\lVert\tilde{\varphi}_{N}\rVert_{h}$ is bounded by

$\lVert\tilde{\varphi}_{N}\rVert_{h}\leqslant\frac{2\left(\frac{hr}{2}\right)^{N}}{N!}\frac{N}{N-1}\leqslant\frac{4\left(\frac{hr}{2}\right)^{N}}{N!}.$

For $N=0$, we have $\lVert\tilde{\varphi}_{N}\rVert_{h}\leqslant
2\lVert\varphi\rVert_{h}\leqslant 2$, which also satisfies (15), thus (15)
holds for all $N\in\mathbb{N}$. ∎

######  Remark 4.

For $\varphi\in\mathcal{S}$, Chebyshev polynomial approximation error
converges with a super-geometric convergence rate. Through orthogonal
projections, it is the best approximation to minimize the norm of the
approximation error [11, 14].

##  3 Main results

In this section, we first present necessary positive semi-definite conditions
stemming from general polynomial approximations of the functional argument
$\varphi$. Second, for the Chebyshev orthogonal polynomials projection, a
sufficiency positive semi-definite condition and a stability condition
verifiable in a tractable number of operations are presented.

###  3.1 Necessary conditions: general polynomial approximations

By substituting the polynomial approximation (10) of $\varphi\in
PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$ into the functional (6), we
obtain the functional $v_{0}(\varphi_{N})$, which is the polynomial
approximation of functional (6). Observe that $I_{1}$ and $I_{2}$ reduce to

$I_{1}(\varphi_{N})+I_{2}(\varphi_{N})=2(\varphi(0)-D\varphi(-h))^{\top}\\\
\times\int_{-h}^{0}\left(U^{\top}(h+\theta)A_{1}-U^{\prime\top}(h+\theta)D\right)\Theta_{N}^{\top}(\theta)\Phi_{N}\mathrm{d}\theta\\\
=2(\varphi(0)-D\varphi(-h))^{\top}(J_{1N}-J_{2N})\Phi_{N},$ (16)

where

$\displaystyle
J_{1N}=\int_{-h}^{0}U^{\top}(h+\theta)A_{1}\Theta_{N}^{\top}(\theta)\mathrm{d}\theta,\quad
J_{2N}=\int_{-h}^{0}U^{\prime\top}(h+\theta)D\Theta_{N}^{\top}(\theta)\mathrm{d}\theta.$

Evaluating the remaining summands of the functional, we get

$\begin{split}\sum\limits_{i=3}^{6}I_{i}(\varphi_{N})&=\int_{-h}^{0}\int_{-h}^{0}\Phi_{N}^{\top}\Theta_{N}(\theta_{1})\left(A_{1}^{\top}U(\theta_{1}-\theta_{2})A_{1}\right.\\\
&\phantom{=}\left.+2A_{1}^{\top}U^{\prime}(\theta_{1}-\theta_{2})D\right)\Theta_{N}^{\top}(\theta_{2})\Phi_{N}\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&\phantom{=}-\int_{-h}^{0}\int_{-h}^{\theta_{1}}\Phi_{N}^{\top}\Theta_{N}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\Theta_{N}^{\top}(\theta_{2})\Phi_{N}\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&\phantom{=}-\int_{-h}^{0}\int_{\theta_{1}}^{0}\Phi_{N}^{\top}\Theta_{N}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\Theta_{N}^{\top}(\theta_{2})\Phi_{N}\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&\phantom{=}-\int_{-h}^{0}\Phi_{N}^{\top}\Theta_{N}(\theta)D^{\top}PD\Theta_{N}^{\top}(\theta)\Phi_{N}\mathrm{d}\theta\\\
&=\Phi_{N}^{\top}(J_{3N}+2J_{4N}-J_{5N}-J_{6N})\Phi_{N}.\end{split}$ (17)

Here,

$\displaystyle J_{3N}$
$\displaystyle=\int_{-h}^{0}\int_{-h}^{0}\Theta_{N}(\theta_{1})A_{1}^{\top}U(\theta_{1}-\theta_{2})A_{1}\Theta_{N}^{\top}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1},$
$\displaystyle J_{4N}$
$\displaystyle=\int_{-h}^{0}\int_{-h}^{0}\Theta_{N}(\theta_{1})A_{1}^{\top}U^{\prime}(\theta_{1}-\theta_{2})D\Theta_{N}^{\top}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1},$
$\displaystyle J_{5N}$
$\displaystyle=\int_{-h}^{0}\left(\int_{-h}^{\theta_{1}}\Theta_{N}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\Theta_{N}^{\top}(\theta_{2})\mathrm{d}\theta_{2}\right.$
$\displaystyle\phantom{=}\left.+\int_{\theta_{1}}^{0}\Theta_{N}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\Theta_{N}^{\top}(\theta_{2})\mathrm{d}\theta_{2}\right)\mathrm{d}\theta_{1},$
$\displaystyle J_{6N}$
$\displaystyle=\int_{-h}^{0}\Theta_{N}(\theta)D^{\top}PD\Theta_{N}^{\top}(\theta)\mathrm{d}\theta.$

Note that, for any $\Psi(\theta)=\Psi^{\top}(-\theta),~{}\theta>0$, we have

$\int\limits_{-h}^{0}\int\limits_{\theta_{1}}^{0}\theta_{1}^{i}\theta_{2}^{j}\Psi(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}=\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\theta_{2}^{i}\theta_{1}^{j}\Psi^{\top}(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}.$

Hence, $J_{3N},~{}J_{4N}$ and$~{}J_{5N}$ can be rewritten as

$\displaystyle J_{3N}$
$\displaystyle=\mathrm{He}\\!\left(\\!\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\Theta_{N}(\theta_{1})A_{1}^{\top}U(\theta_{1}-\theta_{2})A_{1}\Theta_{N}^{\top}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\!\right)\\!,$
$\displaystyle J_{4N}$
$\displaystyle=\mathrm{He}\\!\left(\\!\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\Theta_{N}(\theta_{1})A_{1}^{\top}U^{\prime}(\theta_{1}-\theta_{2})D\Theta_{N}^{\top}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\!\right)\\!,$
$\displaystyle J_{5N}$
$\displaystyle=\mathrm{He}\\!\left(\\!\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\Theta_{N}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\Theta_{N}^{\top}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\!\right)\\!.$

Gathering the approximated summands, we rewrite the polynomial approximation
of functional (6) as:

$v_{0}(\varphi_{N})=\Xi^{\top}\mathbf{P}_{N}\Xi,$ (18)

where

$\mathbf{P}_{N}=\begin{bmatrix}J_{0N}&J_{1N}+J_{2N}\\\
\star&J_{3N}+2J_{4N}-J_{5N}-J_{6N}\end{bmatrix},\quad J_{0N}=U(0),$
$\Xi=[(\varphi(0)-D\varphi(-h))^{\top}\quad\Phi_{N}^{\top}]^{\top}.$

Now, for any polynomial approximation of $\varphi$ and any order of
approximation $N$, we give a necessary stability condition through the
functional approximation (18).

######  Theorem 2.

If system (1) is exponentially stable, then the matrix $\mathbf{P}_{N}\geq
0,~{}\forall N.$

###### Proof.

Since system (1) is exponentially stable, it is also stable, thus by
construction of the functional, we have

$v_{0}(\varphi)=\int_{0}^{+\infty}x^{T}(t,\varphi)Wx(t,\varphi)\textup{d}t\geqslant
0,$

for any $\varphi\in PC^{1}\bigl{(}[-h,0],\mathbb{R}^{n}\bigr{)}$. In
particular, for the polynomial approximation
$\varphi_{N}:[-h,0]\to\mathbb{R}^{n}$ given by (10), it follows from (18) that

$v_{0}(\varphi_{N})=\Xi^{\top}\mathbf{P}_{N}\Xi\geqslant 0,$

for any $\Xi$ and $N$, which implies $\mathbf{P}_{N}\geq 0.$ ∎

######  Remark 5.

Notice that an approximation of any order, even small, provides a necessary
condition.

###  3.2 Sufficient condition: Chebyshev polynomial projection

Here, we focus on the Chebyshev polynomial projection introduced in Section
2.2.3. To conclude on the exponential stability of system (1), we first
analyze how the approximation error of the Chebyshev orthogonal polynomial
affects the functional approximation error $\Upsilon_{N}$. Next, we combine
this result with the sufficient stability condition for functional
$v(\varphi)$ established in Lemma 1.

The approximation error on the functional argument $\varphi$ is conveyed to
functional (6), leading to the following expression for the functional
approximation error

$\Upsilon_{N}=v_{0}(\varphi)-v_{0}(\varphi_{N}),\quad\varphi\in\mathcal{S}.$
(19)

The value $\Upsilon_{N}$ quantifies the error between $v_{0}(\varphi_{N})$ and
$v_{0}(\varphi)$ on the set of functions $\varphi\in\mathcal{S}$. We now
proceed to compute and estimate the functional approximation error in order to
transfer the convergence properties of Chebyshev to the error $\Upsilon_{N}$.

Taking into account that if the Lyapunov condition in Theorem 1 holds, then
$U(s)$ is guaranteed to be bounded, we introduce the constants

$M_{1}=\sup_{\theta\in(0,h)}\|U^{\top}(\theta)A_{1}-U^{\prime\top}(\theta)D\|,$
$M_{2}=\sup_{\theta\in(0,h)}\|A_{1}^{\top}U(\theta)A_{1}+2A_{1}^{\top}U^{\prime}(\theta)D-D^{\top}U^{\prime\prime}(\theta)D\|,$
$M_{3}=\|D^{\top}PD\|.$

Denote
$\Psi(\theta)=A_{1}^{\top}U(\theta)A_{1}+2A_{1}^{\top}U^{\prime}(\theta)D$ and
notice that the term $\Upsilon_{N}$ can be presented as

$\Upsilon_{N}=\sum\limits_{i=1}^{6}I_{i}^{\mathrm{error}},$

where

$I_{1}^{\mathrm{error}}+I_{2}^{\mathrm{error}}=2[\varphi(0)-D\varphi(-h)]^{\top}\int_{-h}^{0}\left(U^{\top}(h+\theta)A_{1}-U^{\prime\top}(h+\theta)D\right)\tilde{\varphi}_{N}(\theta)\mathrm{d}\theta,$
$\begin{split}&\sum\limits_{i=3}^{6}I_{i}^{\mathrm{error}}=\int_{-h}^{0}\int_{-h}^{0}\varphi^{\top}(\theta_{1})\Psi(\theta_{1}-\theta_{2})\varphi(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-\int_{-h}^{0}\int_{-h}^{0}\varphi^{\top}_{N}(\theta_{1})\Psi(\theta_{1}-\theta_{2})\varphi_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-\int_{-h}^{0}\int_{-h}^{\theta_{1}}\varphi^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\varphi(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&+\int_{-h}^{0}\int_{-h}^{\theta_{1}}\varphi_{N}^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\varphi_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-\int_{-h}^{0}\int_{\theta_{1}}^{0}\varphi^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\varphi(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\end{split}$
$\begin{split}&+\int_{-h}^{0}\int_{\theta_{1}}^{0}\varphi_{N}^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\varphi_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-\int_{-h}^{0}\varphi^{\top}(\theta)D^{\top}PD\varphi(\theta)\mathrm{d}\theta+\int_{-h}^{0}\varphi_{N}^{\top}(\theta)D^{\top}PD\varphi_{N}(\theta)\mathrm{d}\theta.\end{split}$

Substituting $\varphi_{N}=\varphi-\tilde{\varphi}_{N}$, we have

$\begin{split}&\sum\limits_{i=3}^{6}I_{i}^{\mathrm{error}}=2\int_{-h}^{0}\int_{-h}^{0}\varphi^{\top}(\theta_{1})\Psi(\theta_{1}-\theta_{2})\tilde{\varphi}_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-\int_{-h}^{0}\int_{-h}^{0}\tilde{\varphi}_{N}^{\top}(\theta_{1})\Psi(\theta_{1}-\theta_{2})\tilde{\varphi}_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-2\int_{-h}^{0}\int_{-h}^{\theta_{1}}\varphi^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\tilde{\varphi}_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&+\int_{-h}^{0}\int_{-h}^{\theta_{1}}\tilde{\varphi}_{N}^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\tilde{\varphi}_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-2\int_{-h}^{0}\int_{\theta_{1}}^{0}\varphi^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\tilde{\varphi}_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&+\int_{-h}^{0}\int_{\theta_{1}}^{0}\tilde{\varphi}_{N}^{\top}(\theta_{1})D^{\top}U^{\prime\prime}(\theta_{1}-\theta_{2})D\tilde{\varphi}_{N}(\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
&-2\int_{-h}^{0}\tilde{\varphi}_{N}^{\top}(\theta)D^{\top}PD\varphi(\theta)\mathrm{d}\theta+\int_{-h}^{0}\tilde{\varphi}_{N}^{\top}(\theta)D^{\top}PD\tilde{\varphi}_{N}(\theta)\mathrm{d}\theta.\end{split}$

Considering $\varphi\in\mathcal{S}$,
$\|\varphi(\theta)\|\leqslant\|\varphi(0)\|=1$ and Lemma 2, the summands of
$\Upsilon_{N}$ are estimated by:

$|I_{1}^{\mathrm{error}}+I_{2}^{\mathrm{error}}|\leqslant
2h(1+\|D\|)M_{1}\|\tilde{\varphi}_{N}\|_{h},$
$\left|\sum\limits_{i=3}^{6}I_{i}^{\mathrm{error}}\right|\leqslant
2h(hM_{2}+M_{3})\|\tilde{\varphi}_{N}\|_{h}+h(hM_{2}+M_{3})\|\tilde{\varphi}_{N}\|_{h}^{2}.$

Hence, the functional approximation error $\Upsilon_{N}$ admits the upper
bound

$|\Upsilon_{N}|\leqslant
2c_{1}\|\tilde{\varphi}_{N}\|_{h}+c_{2}\|\tilde{\varphi}_{N}\|_{h}^{2},$ (20)
$c_{1}=h((1+\|D\|)M_{1}+hM_{2}+M_{3}),\quad c_{2}=h(hM_{2}+M_{3}).$

Using the upper bound (15) on the argument error of the Chebyshev
approximation of $\varphi$, we obtain

$\small|\Upsilon_{N}|\leqslant
2c_{1}\frac{4\left(\frac{hr}{2}\right)^{N}}{N!}+c_{2}\left(\frac{4\left(\frac{hr}{2}\right)^{N}}{N!}\right)^{2}\leqslant\delta_{N}=\frac{(8c_{1}\\!+\\!16c_{2})\mu^{N}}{N!},$
(21)

where $\mu=\max\left\\{\frac{hr}{2},\left(\frac{hr}{2}\right)^{2}\right\\}$
and $r$ is defined in (8).

Next, a sufficiency stability condition based on the functional approximation
(18).

######  Theorem 3.

Let $\mathbf{P}_{N}$ in (18) where $\Phi_{N}$ is determined through Chebyshev
polynomial projection. If the Lyapunov condition holds and there exists a
natural number $N$ such that

$\mathbf{P}_{N}-\mathcal{X}_{N}\geq 0,$ (22)

where the matrix $\mathcal{X}_{N}$ is of the same dimension as
$\mathbf{P}_{N}$ with the upper left $n\times n$ block equal to
$\frac{\delta_{N}}{(1-\|D\|)^{2}}I_{n}$ and all other blocks equal to zero,
and $\delta_{N}$ defined in (21), then system (1) is exponentially stable.

###### Proof.

By contradiction, we assume that system (1) is unstable, but condition (22)
holds. Now, it follows from (18), (19) and (21) that functional (6) admits the
lower bound

$\begin{split}v_{0}(\varphi)&\geqslant\Xi^{\top}\mathbf{P}_{N}\Xi-\delta_{N}\|\varphi(0)\|^{2}\\\
&\geqslant\Xi^{\top}\mathbf{P}_{N}\Xi-\frac{\delta_{N}}{(1-\|D\|)^{2}}\|\varphi(0)-D\varphi(-h)\|^{2}\\\
&=\Xi^{\top}(\mathbf{P}_{N}-\mathcal{X}_{N})\Xi.\end{split}$

According to the assumption of instability of the system and Lemma 1, there
exists a constant $a_{0}>0$ such that

$\Xi^{\top}(\mathbf{P}_{N}-\mathcal{X}_{N})\Xi\leqslant
v_{0}(\varphi)<-a_{0}.$

This implies that $\mathbf{P}_{N}-\mathcal{X}_{N}<0$, which is a
contradiction, thus (22) is proved. ∎

######  Corollary 1.

System (1) is exponentially stable, if and only if the Lyapunov condition
holds and the matrix $\mathbf{P}_{N^{\star}}\geq 0$, where the order of
approximation is given by

$N^{\star}=\left\lceil\mu\exp\left[1+\mathcal{W}\left(\frac{1}{\exp(1)\mu}\log\left(\frac{8c_{1}+16c_{2}}{a_{0}}\right)\right)\right]\right\rceil.$
(23)

Here, the constant $a_{0}$ is determined in Theorem 1.

###### Proof.

[Necessity] The necessity follows from Theorem 2. In particular,
$N=N^{\star}$.

[Sufficiency] By contradiction, assume that system (1) is unstable but that
$\mathbf{P}_{N^{\star}}$ is positive semi-definite. Moreover, recall that
$v_{0}(\varphi)$ can be expressed as

$v_{0}(\varphi)=v_{0}(\varphi_{N^{\star}})+\Upsilon_{N^{\star}}.$

By Lemma 1, there exists $\varphi\in\mathcal{S}$ such that
$v_{0}(\varphi)<-a_{0}$, which implies that

$v_{0}(\varphi_{N^{\star}})=v_{0}(\varphi)-\Upsilon_{N^{\star}}<-a_{0}+|\Upsilon_{N^{\star}}|.$

Then, for $N^{\star}$, it follows from (20) that the functional approximation
error admits the upper bound

$|\Upsilon_{N^{\star}}|\leqslant\frac{(8c_{1}+16c_{2})\mu^{N^{\star}}}{N^{\star}!}.$

Applying the logarithm to (20) and the Maclaurin integral test give the
following upper bound

$\begin{split}\log\left(|\Upsilon_{N^{\star}}|\right)&\leqslant\log(8c_{1}+16c_{2})+N^{\star}\log\left(\mu\right)-\int_{0}^{N^{\star}}\log(s)\mathrm{d}s\\\
&=\log\left(8c_{1}+16c_{2}\right)-N^{\star}\log\left(\frac{N^{\star}}{\exp(1)\mu}\right),\end{split}$

equivalently,

$\frac{1}{\exp(1)\mu}\log\left(\frac{|\Upsilon_{N^{\star}}|}{8c_{1}+16c_{2}}\right)\leqslant-\frac{N^{\star}}{\exp(1)\mu}\log\left(\frac{N^{\star}}{\exp(1)\mu}\right).$

Under the condition

$\log\left(\frac{N}{\exp(1)\mu}\right)\geqslant\mathcal{W}\left(\frac{1}{\exp(1)\mu}\log\left(\frac{8c_{1}+16c_{2}}{a_{0}}\right)\right),$

which is satisfied for order $N^{\star}$ given by (23), we have
$\log(|\Upsilon_{N^{\star}}|)\leqslant\log(a_{0})$ implying that
$|\Upsilon_{N^{\star}}|\leqslant a_{0}$. Hence,

$v_{0}(\varphi_{N^{\star}})<-a_{0}+|\Upsilon_{N^{\star}}|\leqslant 0,$

contradicting the statement of positive semi-definiteness. ∎

Corollary 1 merges Theorem 2 and 3 to reach the following necessary and
sufficient stability condition.

##  4 Recursive method for projection computation

Following closely the ideas introduced in [6], we devote this section to
present efficient tools for computing the integrals involved in determining
matrix $\mathbf{P}_{N}$. These matrices are essential for conducting the
stability test discussed in the preceding sections. The computation of these
integrals imposes a significant numerical burden, potentially compromising the
proper evaluation of the condition. Our aim is to ensure accurate condition
testing by reducing the computational load.

Observe that $J_{1N},~{}J_{2N},~{}J_{3N},~{}J_{4N},~{}J_{5N}$ and $J_{6N}$
defined in (16) and (17) involve the following elementary integral matrices
for $k=\overline{0,N-1}$ and $i,j=\overline{0,N-1}$:

$\displaystyle
G_{k}=\int\limits_{-h}^{0}U(h+\theta)\theta^{k}\mathrm{d}s,\quad\bar{G}_{k}=\int\limits_{-h}^{0}U(\theta)\theta^{k}\mathrm{d}\theta,$
$\displaystyle
H_{ij}=\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1},\quad\bar{H}_{ij}=\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2}-h)\mathrm{d}\theta_{2}\mathrm{d}\theta_{1},$

which must be computed. To carry out this task, let us consider the dynamic
property (2), denote

$L=\begin{bmatrix}I\otimes I&-I\otimes D\\\ -D^{\top}\otimes I&I\otimes
I\end{bmatrix}^{-1}\begin{bmatrix}I\otimes A_{0}&I\otimes A_{1}\\\
-A_{1}^{\top}\otimes I&-A_{0}^{\top}\otimes I\end{bmatrix},$

and define the vectors $G_{k}^{v}=\textup{vec}(G_{k})$,
$\bar{G}_{k}^{v}=\textup{vec}(\bar{G}_{k})$,
$H_{ij}^{v}=\textup{vec}(H_{ij})$,
$\bar{H}_{ij}^{v}=\textup{vec}(\bar{H}_{ij})$, $U^{v}(0)=\textup{vec}(U(0))$
and $U^{v}(-h)=\textup{vec}(U(-h))$. Next, we introduce two propositions
allowing the recursive computation of the above integral.

######  Proposition 1.

If the Lyapunov condition holds and $\mathrm{det}(L)\neq 0$, then the set of
matrices $\\{G_{k}\\}_{k=\overline{0,N-1}}$ can be computed as

$\displaystyle\begin{bmatrix}G_{0}^{v}\\\ \bar{G}_{0}^{v}\end{bmatrix}$
$\displaystyle=L^{-1}\begin{bmatrix}U^{v}(h)-U^{v}(0)\\\
U^{v}(0)-U^{v}(-h)\end{bmatrix},$ (24)
$\displaystyle\begin{bmatrix}G_{k}^{v}\\\ \bar{G}_{k}^{v}\end{bmatrix}$
$\displaystyle=-L^{-1}\left(k\begin{bmatrix}G_{k-1}^{v}\\\
\bar{G}_{k-1}^{v}\end{bmatrix}+(-h)^{k}\begin{bmatrix}U^{v}(0)\\\
U^{v}(-h)\end{bmatrix}\right),\,k\neq 0.$ (25)

###### Proof.

Integrating the dynamic property (2) yields

$\int\limits_{-h}^{0}U^{\prime}(h+\theta)\mathrm{d}\theta-\int\limits_{-h}^{0}U^{\prime}(\theta)\mathrm{d}\theta
D=\int\limits_{-h}^{0}U(h+\theta)\mathrm{d}\theta
A_{0}+\int\limits_{-h}^{0}U(\theta)\mathrm{d}\theta A_{1}.$

We also note that

$\int_{-h}^{0}U(h+\theta)\mathrm{d}\theta=G_{0},\quad\bar{G}_{0}=\int\limits_{-h}^{0}U(\theta)\mathrm{d}\theta=\int\limits_{-h}^{0}U^{\top}(-\theta)\mathrm{d}\theta=\int\limits_{0}^{h}U^{\top}(\theta)\mathrm{d}\theta=\int\limits_{-h}^{0}U^{\top}(h+\theta)\mathrm{d}\theta=G_{0}^{\top}.$

Hence, we have

$U(h)-U(0)-U(0)D+U(-h)D=G_{0}A_{0}+G_{0}^{\top}A_{1}.$

Notice that this equation can be rewritten as

$U(0)-U(-h)-D^{\top}U(h)D+D^{\top}U(0)=-A_{0}^{\top}G_{0}^{\top}-A_{1}^{\top}G_{0}.$

By vectorization based on Kronecker product properties [17], we arrive at the
following system of linear algebraic equations

$\begin{bmatrix}I\otimes A_{0}&I\otimes A_{1}\\\ -A_{1}^{\top}\otimes
I&-A_{0}^{\top}\otimes I\end{bmatrix}\begin{bmatrix}G_{0}^{v}\\\
\bar{G}_{0}^{v}\end{bmatrix}=\begin{bmatrix}I\otimes I&-I\otimes D\\\
-D^{\top}\otimes I&I\otimes I\end{bmatrix}\begin{bmatrix}U^{v}(h)-U^{v}(0)\\\
U^{v}(0)-U^{v}(-h)\end{bmatrix}.$

Since $\mathrm{det}(L)\neq 0$, we arrive at a unique solution of this system
given by (24).
Following the same idea, the dynamic property (2) is multiplied by
$\theta^{k}$ and integrated to obtain

$\int\limits_{-h}^{0}U^{\prime}(h+\theta)\theta^{k}\mathrm{d}\theta-\int\limits_{-h}^{0}U^{\prime}(\theta)\theta^{k}\mathrm{d}\theta
D=\int\limits_{-h}^{0}U(h+\theta)\theta^{k}\mathrm{d}\theta
A_{0}+\int\limits_{-h}^{0}U(\theta)\theta^{k}\mathrm{d}\theta A_{1}.$

Applying integration by parts leads to

$-kG_{k-1}-(-h)^{k}U(0)+k\bar{G}_{k-1}D+(-h)^{k}U(-h)D=G_{k}A_{0}+\bar{G}_{k}A_{1}.$

Similar steps to the dynamic property (5) give the two following algebraic
equations rewritten in vectorized form

$\begin{bmatrix}I\otimes A_{0}&I\otimes A_{1}\\\ -A_{1}^{\top}\otimes
I&-A_{0}^{\top}\otimes I\end{bmatrix}\begin{bmatrix}G_{k}^{v}\\\
\bar{G}_{k}^{v}\end{bmatrix}=-\begin{bmatrix}I\otimes I&-I\otimes D\\\
-D^{\top}\otimes I&I\otimes I\end{bmatrix}\left(k\begin{bmatrix}G_{k-1}^{v}\\\
\bar{G}_{k-1}^{v}\end{bmatrix}+(-h)^{k}\begin{bmatrix}U^{v}(0)\\\
U^{v}(-h)\end{bmatrix}\right).$

equivalently,

$L\begin{bmatrix}G_{k}^{v}\\\
\bar{G}_{k}^{v}\end{bmatrix}=-\left(k\begin{bmatrix}G_{k-1}^{v}\\\
\bar{G}_{k-1}^{v}\end{bmatrix}+(-h)^{k}\begin{bmatrix}U^{v}(0)\\\
U^{v}(-h)\end{bmatrix}\right).$

When $\mathrm{det}(L)\neq 0$, we conclude that the unique solution is given by
(25). ∎

######  Proposition 2.

If the Lyapunov condition holds and $\mathrm{det}(L)\neq 0$, then the set of
matrices $\\{H_{ij}\\}_{i,j=\overline{0,N-1}}$ can be computed as

$\begin{bmatrix}H_{ij}^{v}\\\
\bar{H}_{ij}^{v}\end{bmatrix}=L^{-1}\begin{pmatrix}\frac{(-h)^{i+j+1}}{i+j+1}U^{v}(0)+(-h)^{j}G_{i}^{v}+jH_{i,j-1}^{v}\\\
\frac{(-h)^{i+j+1}}{i+j+1}U^{v}(-h)+(-h)^{j}\bar{G}_{i}^{v}+j\bar{H}_{i,j-1}^{v}\end{pmatrix}.$
(26)

###### Proof.

Repeating the steps used in the proof of Proposition 1, the dynamic property
(2) for $\theta_{1}-\theta_{2}>0$ gives

$-\frac{\partial U(\theta_{1}-\theta_{2})}{\partial\theta_{2}}+\frac{\partial
U(\theta_{1}-\theta_{2}-h)}{\partial\theta_{2}}D=U(\theta_{1}-\theta_{2})A_{0}+U(\theta_{1}-\theta_{2}-h)A_{1}.$

Now, pre and post multiply the discretized dynamic property (2) by
$\theta_{1}^{i}$ and $\theta_{2}^{j}$, respectively:

$-\theta_{1}^{i}\theta_{2}^{j}\frac{\partial
U(\theta_{1}-\theta_{2})}{\partial\theta_{2}}+\theta_{1}^{i}\theta_{2}^{j}\frac{\partial
U(\theta_{1}-\theta_{2}-h)}{\partial\theta_{2}}D=\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2})A_{0}+\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2}-h)A_{1}.$

Then, we integrate with respect to $\theta_{2}$ the previous expression, from
$-h$ up to $\theta_{1}$ in such a way that the difference
$\theta_{1}-\theta_{2}>0$. It implies that

$-\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}\frac{\partial
U(\theta_{1}-\theta_{2})}{\partial\theta_{2}}\mathrm{d}\theta_{2}+\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}\frac{\partial
U(\theta_{1}-\theta_{2}-h)}{\partial\theta_{2}}\mathrm{d}\theta_{2}D=\\\
\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}A_{0}+\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2}-h)\mathrm{d}\theta_{2}A_{1}.$

Compute the integrals in the above expression:

$\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}\frac{\partial
U(\theta_{1}-\theta_{2})}{\partial\theta_{2}}\mathrm{d}\theta_{2}=\theta_{1}^{i+j}U(0)-(-h)^{j}\theta_{1}^{i}U(h+\theta_{1})-\int\limits_{-h}^{\theta_{1}}j\theta_{1}^{i}\theta_{2}^{j-1}U(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2},$
$\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}\frac{\partial
U(\theta_{1}-\theta_{2}-h)}{\partial\theta_{2}}\mathrm{d}\theta_{2}=\theta_{1}^{i+j}U^{\top}(h)-(-h)^{j}\theta_{1}^{i}U(\theta_{1})-\int\limits_{-h}^{\theta_{1}}j\theta_{1}^{i}\theta_{2}^{j-1}U(\theta_{1}-\theta_{2}-h)\mathrm{d}\theta_{2}.$

noting that the last term vanishes in the case $j=0$. Thus,

$-\theta_{1}^{i+j}U(0)+(-h)^{j}\theta_{1}^{i}U(h+\theta_{1})+\int\limits_{-h}^{\theta_{1}}j\theta_{1}^{i}\theta_{2}^{j-1}U(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}+\theta_{1}^{i+j}U^{\top}(h)D-(-h)^{j}\theta_{1}^{i}U(\theta_{1})D\\\
-\int\limits_{-h}^{\theta_{1}}j\theta_{1}^{i}\theta_{2}^{j-1}U(\theta_{1}-\theta_{2}-h)\mathrm{d}\theta_{2}D=\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}A_{0}+\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2}-h)\mathrm{d}\theta_{2}A_{1}.$

Further, integrate with respect to $\theta_{1}$

$-\int\limits_{-h}^{0}\theta_{1}^{i+j}U(0)\mathrm{d}\theta_{1}+\int\limits_{-h}^{0}(-h)^{j}\theta_{1}^{i}U(h+\theta_{1})\mathrm{d}\theta_{1}+\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}j\theta_{1}^{i}\theta_{2}^{j-1}U(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}\\\
+\int\limits_{-h}^{0}\theta_{1}^{i+j}U^{\top}(h)D\mathrm{d}\theta_{1}-\int\limits_{-h}^{0}(-h)^{j}\theta_{1}^{i}U(\theta_{1})\mathrm{d}\theta_{1}D-\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}j\theta_{1}^{i}\theta_{2}^{j-1}U(\theta_{1}-\theta_{2}-h)\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}D=\\\
\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2})\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}A_{0}+\int\limits_{-h}^{0}\int\limits_{-h}^{\theta_{1}}\theta_{1}^{i}\theta_{2}^{j}U(\theta_{1}-\theta_{2}-h)\mathrm{d}\theta_{2}\mathrm{d}\theta_{1}A_{1},$

equivalently,

$\frac{(-h)^{i+j+1}}{i+j+1}(U(0)-U^{\top}(h)D)+(-h)^{j}G_{i}+jH_{i,j-1}-(-h)^{j}\bar{G}_{i}D-j\bar{H}_{i,j-1}D=H_{ij}A_{0}+\bar{H}_{ij}A_{1}.$

Using similar steps for the dynamic property (5), we arrive at the two
following algebraic equations

$\frac{(-h)^{i+j+1}}{i+j+1}(U(0)-U^{\top}(h)D)+(-h)^{j}G_{i}+jH_{i,j-1}-(-h)^{j}\bar{G}_{i}D-j\bar{H}_{i,j-1}D=H_{ij}A_{0}+\bar{H}_{ij}A_{1},$
$\frac{(-h)^{i+j+1}}{i+j+1}(U^{\top}(h)-D^{\top}U(0))+(-h)^{j}\bar{G}_{i}+j\bar{H}_{i,j-1}-D^{\top}(-h)^{j}\bar{G}_{i}-D^{\top}jH_{i,j-1}=-A_{1}^{\top}H_{ij}-A_{0}^{\top}\bar{H}_{ij}.$

By vectorization, the following system of linear algebraic equations is
obtained

$\begin{bmatrix}I\otimes A_{0}&I\otimes A_{1}\\\ -A_{1}^{\top}\otimes
I&-A_{0}^{\top}\otimes I\end{bmatrix}\begin{bmatrix}H_{ij}^{v}\\\
\bar{H}_{ij}^{v}\end{bmatrix}=\begin{bmatrix}I\otimes I&-I\otimes D\\\
-D^{\top}\otimes I&I\otimes
I\end{bmatrix}\begin{bmatrix}\frac{(-h)^{i+j+1}}{i+j+1}U(0)+(-h)^{j}G_{i}^{v}+jH_{i,j-1}^{v}\\\
\frac{(-h)^{i+j+1}}{i+j+1}U^{v}(-h)+(-h)^{j}\bar{G}_{i}^{v}+j\bar{H}_{i,j-1}^{v}\end{bmatrix}.$

When $\mathrm{det}(L)\neq 0$, we conclude that the unique solution is given by
(26). ∎

######  Remark 6.

Matrices $(G_{k},\bar{G}_{k},H_{ij},\bar{H}_{ij})$ are finally obtained by
devectorization of $(G_{k}^{v},\bar{G}_{k}^{v},H_{ij}^{v},\bar{H}_{ij}^{v})$.

With the help of Proposition 1 and 2, we now give an explicit expression for
the integral entries of matrix $\mathbf{P}_{N}$. Thus, the integrals involving
the delay Lyapunov matrix and its derivatives are given by

$\displaystyle J_{0N}$ $\displaystyle=U(0),$ $\displaystyle J_{1N}$
$\displaystyle=\begin{bmatrix}G_{0}^{\top}A_{1}&\cdots&G_{N-1}^{\top}A_{1}\end{bmatrix},$
$\displaystyle J_{2N}$
$\displaystyle=-\begin{bmatrix}-U^{\top}(h)&G_{0}^{\top}&\\!\cdots\\!&(N-1)G_{N-2}^{\top}\end{bmatrix}D-\begin{bmatrix}1&\\!\cdots\\!&(-h)^{N-1}\end{bmatrix}U(0)D,$
$\displaystyle J_{3N}$
$\displaystyle=\mathrm{He}\left(\left\\{A_{1}^{\top}H_{ij}A_{1}\right\\}_{i,j=0}^{N-1}\right),$
$\displaystyle J_{4N}$
$\displaystyle=\mathrm{He}\left(\left\\{A_{1}^{\top}\begin{pmatrix}\frac{(-h)^{i+j+1}}{i+j+1}U(0)+(-h)^{j}G_{i}+jH_{i,j-1}\end{pmatrix}D\right\\}_{i,j=0}^{N-1}\right),$
$\displaystyle J_{5N}\\!+\\!J_{6N}$
$\displaystyle=\mathrm{He}\left(\left\\{D^{\top}\begin{pmatrix}-(-h)^{i+j}U(0)-i(-h)^{j}G_{i-1}\end{pmatrix}D\right\\}_{i,j=0}^{N-1}\right)$
$\displaystyle+\mathrm{He}\left(\left\\{D^{\top}j\begin{pmatrix}\frac{(-h)^{i+j}}{i+j}U(0)+(-h)^{j-1}G_{i}+(j-1)H_{i,j-2}\end{pmatrix}D\right\\}_{i,j=0}^{N-1}\right).$

##  5 Illustrative examples

Two examples validate the condition of Corollary 1. The delay Lyapunov matrix
involved in each criterion is computed via the semi-analytic method introduced
in [2] for $W=I_{n}$. The positive semi-definiteness of $\mathbf{P}_{N}$ is
verified by using the function “cholcov” in Matlab as suggested in [6]. In
each example, the stability boundaries obtained by the D-subdivision are
depicted by a solid line. For comparison purposes, we present the order of
approximations $\hat{r}$ and $N_{0}$ introduced in [4] and [5], respectively.

Example 1: Consider the scalar neutral type system

$\frac{\mathrm{d}}{\mathrm{d}t}[x(t)-dx(t-h)]=a_{0}x(t)+a_{1}x(t-h).$

The map of tractable approximation orders for the stable region of the space
of parameter $(a_{0},a_{1})$, shown in Figure 1 where the color code is
related to the order $N^{\star}$ of $\mathbf{P}_{N^{\star}}$, illustrates the
stability condition of Corollary 1.

Figure 1: Map of the order $N^{\star}$ in the space of parameter
$(a_{0},a_{1})$.

For implementing the recursive method described in Section 4, the bits of
precision in Matlab play an important role. To show this, the non-negativity
of $\mathbf{P}_{N}$ is verified by computing its minimum eigenvalue. In Table
1, the minimum eigenvalue is computed for $8,~{}16$ and $32$ bits and
parameter values $(a_{0},a_{1})=(0.8,-1.2)$, $d=-0.3$ and dimension $N=20$. It
is worth mentioning that, see Figure 1, system (5) is exponentially stable for
these parameters. Notice that the validation of the non-negativity of
$\mathbf{P}_{N}$ fails for $8$ and $16$ bits. Thus, the bits of precision must
be increased to implement the recursive method, impacting the execution time
for the validation of the stability condition $\mathbf{P}_{N}$.

Table 1: Minimum eigenvalue of $\mathbf{P}_{20}$ with respect to the precision Bits of precision |  $8$ bits |  $16$ bits |  $32$ bits
---|---|---|---
$\lambda_{\min}(\mathbf{P}_{20})$ | $-4.8^{23}$ | $-2.8^{-19}$ | $1.6^{-26}$

Example 2: The $\sigma$-stability analysis of the proportional-integral
control of a passive linear system leads to studying a quasipolynomial of
neutral type [18]. Its time domain representation is of the form (1), with
matrices $D=\begin{pmatrix}0&0\\\
0&-\frac{\alpha_{2}}{\alpha_{1}}\end{pmatrix}$,

$A_{0}=\frac{1}{\alpha_{1}}\begin{pmatrix}0&\alpha_{1}\\\
-\sigma^{2}\alpha_{1}+\sigma\beta_{1}-\gamma_{1}&-\beta_{1}+2\sigma\alpha_{1}\end{pmatrix},$

$A_{1}=\frac{1}{\alpha_{1}}\begin{pmatrix}0&0\\\
-\sigma^{2}\alpha_{2}+\sigma\beta_{2}-\gamma_{2}&-\beta_{2}+2\sigma\alpha_{2}\end{pmatrix},$

where

$\begin{split}\alpha_{1}&=d+k_{p},\quad\gamma_{1}=bk_{i}d+ak_{i},\\\
\alpha_{2}&=(d-k_{p})e^{\sigma h},\quad\gamma_{2}=(bk_{i}d-ak_{i})e^{\sigma
h},\\\ \beta_{1}&=(bk_{p}+a)d+bd^{2}+ak_{p}+k_{i},\\\
\beta_{2}&=((bk_{p}+a)d-bd^{2}-ak_{p}-k_{i})e^{\sigma h}.\end{split}$

The stability of the difference operator imposes in the D-subdivision map the
additional condition $|k_{p}|<26.67$. For parameter values
$a=0.4,~{}b=50,~{}h=0.2,~{}d=0.8$ and $\sigma=0.3$, the map of the
approximation order $N^{\star}$ for the space of parameter $(k_{p},k_{i})$ is
depicted in Figure 2.

Figure 2: Map of the order $N^{\star}$ in the space of parameter
$(k_{p},k_{i})$.

The order of approximation $N^{\star}$ in Corollary 1 is computed and compared
with the results obtained in [4] and [5] for two different points of the space
of parameters. Table 2 shows that the value $N^{\star}$ is significantly
smaller compared to $N_{0}$ and $\hat{r}$, where $(\textbf{---})$ indicates
that the result is not verifiable computationally due to it exceeds the
computer RAM. The favorable outcome of Corollary 1 stems from using the
Chebyshev polynomials instead of functions based on the fundamental matrix [4]
and piece-wise linear approximation [5], resulting in a faster functional
argument approximation convergence rate. Thus, we can conclude the stability
of neutral-type systems in a small order of approximations.

Table 2: Stability criteria given by Corollary 1, Theorem 6 in [5] and Theorem 4 in [4] $(k_{p},k_{i})$ | $N^{\star}$ | Result |  $N_{0}$ [5] | Result |  $\hat{r}$ [4] | Result
---|---|---|---|---|---|---
$(1,1)$ | $26$ | Stable | $65$ | Stable | $5\times 10^{22}$ | —
$(1,-1)$ | $31$ | Unstable | $89$ | Unstable | $3\times 10^{26}$ | —

##  6 Conclusions

We present necessary stability conditions for neutral time-delay systems in
terms of a positive semi-definite test based on any polynomial approximation
of the functional argument. The condition takes a quadratic form independent
of the approximation coefficients. Its dimension depends on the polynomial
approximation method. The sufficiency of the result is proven for the
particular case of Chebyshev polynomial approximation. An estimation of a
dimension for which the positive semi-definiteness test holds is given. Two
examples show the benefits of using Chebyshev in terms of computational
complexity.

## CRediT authorship contribution statement

Gerson Portilla: Methodology, Software, Validation, Writing - Review &
Editing. Mathieu Bajodek: Conceptualization, Methodology, Writing - Review &
Editing preparation. Sabine Mondié: Conceptualization, Supervision, Writing -
Review & Editing.

## References

  * [1] G. Temple, Stability of Motion. Applications of Lyapunov’s second method to differential systems and equations with delay. By N N. Krasovskii. Translated by J L. Brenner. Standford University Press, 1963. Pp. 1–188. 48s., The Mathematical Gazette 49 (367) (1965) 114–114.  doi:10.1017/S0025557200073836.
  * [2] V. L. Kharitonov, Time-delay systems: Lyapunov functionals and matrices, Birkhäuser, 2013.
  * [3] I. V. Alexandrova, A. P. Zhabko, Stability of neutral type delay systems: A joint Lyapunov–Krasovskii and Razumikhin approach, Automatica 106 (2019) 83–90.
  * [4] M. A. Gomez, A. V. Egorov, S. Mondié, Necessary and sufficient stability condition by finite number of mathematical operations for time-delay systems of neutral type, IEEE Transactions on Automatic Control 66 (6) (2021) 2802–2808.  doi:10.1109/TAC.2020.3008392.
  * [5] G. Portilla, I. V. Alexandrova, S. Mondié, Stability tests for neutral-type delay systems: A delay Lyapunov matrix and piecewise linear approximation of the functional argument approach, International Journal of Robust and Nonlinear Control 34 (9) (2024) 5664–5685.  doi:https://doi.org/10.1002/rnc.7293.
  * [6] I. V. Alexandrova, A finite Lyapunov matrix-based stability criterion for linear delay systems via piecewise linear approximation, Systems & Control Letters 181 (2023) 105624.
  * [7] M. Bajodek, F. Gouaisbaut, A. Seuret, Necessary and sufficient stability condition for time-delay systems arising from Legendre approximation, IEEE Transactions on Automatic Control (2022) 1–7doi:10.1109/TAC.2022.3232052.
  * [8] E. Fridman, Introduction to time-delay systems: Analysis and control, Springer, 2014\.
  * [9] A. Seuret, F. Gouaisbaut, Hierarchy of LMI conditions for the stability analysis of time-delay systems, Systems & Control Letters 81 (2015) 1–7.
  * [10] G. M. Phillips, Interpolation and approximation by polynomials, Vol. 14, Springer Science & Business Media, 2003.
  * [11] J. P. Boyd, Chebyshev and Fourier spectral methods, Courier Corporation, 2001.
  * [12] C. A. Fletcher, Computational Galerkin methods, Springer, 1984.
  * [13] T. H. Scholl, V. Hagenmeyer, L. Gröll, What ODE-approximation schemes of time-delay systems reveal about Lyapunov-Krasovskii functionals, Submitted to Systems & Control Letters.
  * [14] L. Fox, I. B. Parker, Chebyshev polynomials in numerical analysis, Oxford University Press.
  * [15] H. Wang, S. Xiang, On the convergence rates of Legendre approximation, Mathematics of computation 81 (278) (2012) 861–877.
  * [16] L. N. Trefethen, Is Gauss quadrature better than Clenshaw–Curtis?, SIAM review 50 (1) (2008) 67–87.
  * [17] H. Roger, R. J. Charles, Topics in matrix analysis, Cambridge University Press Cambridge, UK, 1994.
  * [18] F. Castaños, E. Estrada, S. Mondié, A. Ramírez, Passivity-based PI control of first-order systems with I/O communication delays: a frequency domain analysis, International Journal of Control 91 (11) (2018) 2549–2562.