# The impact of void-finding algorithms on galaxy classification

Fatima Zaidouni Now at Department of Physics, Massachusetts Institute of
Technology, Cambridge, MA 02139  Department of Physics and Astronomy,
University of Rochester, Rochester, NY USA 14627  Dahlia Veyrat Department of
Physics and Astronomy, University of Rochester, Rochester, NY USA 14627  Kelly
A. Douglass Department of Physics and Astronomy, University of Rochester,
Rochester, NY USA 14627 Kelly A. Douglass<EMAIL_ADDRESS>Segev
BenZvi Department of Physics and Astronomy, University of Rochester,
Rochester, NY USA 14627

###### Abstract

We explore how the definition of a void influences the conclusions drawn about
the impact of the void environment on galactic properties using two void-
finding algorithms in the Void Analysis Software Toolkit: V2, a Python
implementation of ZOBOV, and VoidFinder, an algorithm which grows and merges
spherical void regions. Using the Sloan Digital Sky Survey Data Release 7, we
find that galaxies found in VoidFinder voids tend to be bluer, fainter, and
have higher (specific) star formation rates than galaxies in denser regions.
Conversely, galaxies found in V2 voids show no significant differences when
compared to galaxies in denser regions, inconsistent with the large-scale
environmental effects on galaxy properties expected from both simulations and
previous observations. These results align with previous simulation results
that show V2-identified voids “leaking” into the dense walls between voids
because their boundaries extend up to the density maxima in the walls. As a
result, when using ZOBOV-based void finders, galaxies likely to be part of
wall regions are instead classified as void galaxies, a misclassification that
can be critical to our understanding of galaxy evolution.

##  1 Introduction

The large-scale structure of the universe observed in galaxy redshift surveys
exhibits a cosmic web-like structure (Bond et al., 1996) comprising voids, or
vast underdense regions tens of megaparsecs across, that are bordered by
filaments of galaxies that flow into large galaxy clusters (de Lapparent et
al., 1986). Because of the low matter density within voids, their
gravitational evolution remains in the linear regime (Sutter et al., 2012a),
and the few galaxies found in voids evolve relatively free of interactions
with other galaxies.

Voids provide unique laboratories for studying both cosmology and galaxy
physics. Because their evolution is partly shaped by the expansion of the
universe (e.g., Lavaux & Wandelt, 2012; Pisani et al., 2015; Sahlén, 2019),
voids are also useful for studying dark energy and provide precision
constraints on cosmological models via the Alcock-Paczyński test (e.g., Lavaux
& Wandelt, 2012; Sutter et al., 2012b, 2014; Hamaus et al., 2016; Mao et al.,
2017; Nadathur et al., 2019; Hamaus et al., 2020), as well as in neutrino
studies and models of galaxy formation (Peebles, 2001; Constantin et al.,
2008; Sahlén, 2019). In addition, voids are used in measurements of baryon
acoustic oscillations via void-galaxy correlation functions (Nadathur et al.,
2019; Zhao et al., 2020, 2022). Due to the extreme low density of void
interiors, void galaxies rarely interact with other objects, in contrast to
galaxies found in denser regions. As a result, they provide relatively
undisturbed environments for studying galactic evolution.

While voids do have some universal properties, such as an underdense center
($\lesssim 5\%$ of the mean galaxy density) and a steep density gradient at
their boundary (Colberg et al., 2008), the definition of what is a void is
rather vague. In part, this is because void morphologies and densities exhibit
significant variability. Due to the ambiguity, a number of different void
finding algorithms exist to identify underdense regions.

When making conclusions based on different void catalogs, care must be taken
to untangle physical effects from the details of how voids are defined. The
Aspen-Amsterdam void finder comparison project (Colberg et al., 2008)
performed a preliminary comparison between 13 void-finding algorithms,
examining the distributions of the different algorithms’ void properties. A
recent comparison of the void-finding algorithms in the Void Analysis Software
Toolkit (VAST; Douglass et al., 2022) by Veyrat et al. (2023) on their ability
to recover regions of the large-scale structure whose dynamics match those
expected for voids found that the different algorithms do identify different
dynamical regions of the universe. However, the impact of the choice of void-
finding algorithms on either cosmology or galaxy science has not been
systematically studied. Muldrew et al. (2012) used mock galaxy catalogs to
study how various definitions of the large-scale environment change the
inferred effect of the environment on galaxy evolution. In this work, we focus
on the impact of the void-finding algorithms on inferred galaxy properties by
comparing the properties of the galaxies inside voids defined by different
algorithms, but it does not focus on void definitions or void algorithms. This
study allows us to understand how conclusions about the effects of the void
environment depend on the choice of void-finding algorithm.

##  2 Expected Galaxy Evolution in Voids

N-body simulations and semi-analytic models of $\Lambda$CDM cosmology predict
that void galaxies are retarded in their star formation relative to galaxies
in denser regions (Cen, 2011). Numerous studies have been conducted to test
this with observations; investigating the influence of the void environment on
the galactic mass function (e.g., Hoyle et al., 2005; Croton et al., 2005;
Moorman et al., 2015), on galaxy photometry (e.g., Grogin & Geller, 2000;
Rojas et al., 2004, 2005; Patiri et al., 2006; Park et al., 2007; von Benda-
Beckmann & Müller, 2008; Hoyle et al., 2012), on the star formation rate
(e.g., Rojas et al., 2005; Moorman et al., 2016; Beygu et al., 2016), and on
chemical composition (e.g., Kreckel et al., 2012; Douglass & Vogeley, 2017).
The vast majority of these observational studies confirm that galaxies in
voids are fainter and bluer than galaxies in denser regions. This is likely a
result of the retention of gas by void galaxies as they evolve in an
environment with fewer interactions such as mergers, tidal stripping, ram-
pressure stripping, etc., compared to galaxies in denser regions (Cen, 2011).

We therefore expect void galaxies to be fainter, bluer, to have lower mass,
and to have higher (specific) star formation rates regardless of the void-
finding method used in the analysis. We test this by comparing the properties
of void galaxies as classified using two popular void-finding algorithms. To
achieve high statistical significance, we use a large sample of galaxies from
the Sloan Digital Sky Survey Data Release 7 (SDSS DR7, Abazajian et al.,
2009), which to our knowledge is the largest galaxy sample used to study
properties of void galaxies to date. The comparisons shed light on the impact
of different void-finding algorithms in galactic classification.

##  3 Galaxy Catalog: SDSS DR7

SDSS DR7 (Abazajian et al., 2009) is the final public data release from SDSS-
II and was designed to collect spectra of a magnitude-complete galaxy sample,
observing everything brighter than a Petrosian $r$-band magnitude
$m_{r}<17.77$ (Strauss et al., 2002). This final catalog includes 11,663 deg2
of imaging data and five-band photometry for 357 million objects (Lupton et
al., 2001). It also has complete spectroscopy with a resolution
$\lambda/\Delta\lambda\sim 1800$ in the observed wavelength range 3800Å–9200Å
(Blanton et al., 2003) covering a contiguous 9380 deg2 of the Northern
Galactic Cap. In total, the SDSS DR7 spectroscopic sample contains 1.6 million
spectra including 930,000 galaxies, 120,000 quasars, and 460,000 stars.

We use SDSS DR7 as reprocessed and presented in the NASA-Sloan
Atlas111http://www.nsatlas.org/ (NSA, version 1.0.1), a catalog of images and
parameters derived from the SDSS imaging (Blanton et al., 2011). The NSA
provides $K$-corrected absolute magnitudes, colors, and stellar masses. Star
formation rates and specific star formation rates are taken from the MPA-JHU
value-added galaxy catalog222https://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/
based on the method described in Brinchmann et al. (2004). We assume a flat
$\Lambda$CDM cosmology with $\Omega_{M}=0.315$ and $H_{0}=100h$ km s-1 Mpc-1.

##  4 The Void Analysis Software Toolkit

Figure 1: Results from VAST’s VoidFinder (left) and V2 (right) for a
10$h^{-1}$ Mpc-thick slice centered on declination $\delta=32^{\circ}$ of the
main SDSS DR7 galaxy catalog. The intersection of the midplane of the
declination range and the voids are shown in blue; galaxies that fall within a
void are colored in red, while those that do not are shown in black.
Projection effects can produce void galaxies appearing to fall outside a void,
and some wall galaxies appearing to fall within a void.

We use the void catalogs presented in Douglass et al. (2023). The catalogs
were produced using a volume-limited sample of SDSS DR7 ($M_{r}\leq-20$,
$z\leq 0.114$) with two popular void-finding algorithms implemented in the
Void Analysis Software Toolkit (VAST; Douglass et al., 2022): VoidFinder and
V2. VoidFinder (El-Ad & Piran, 1997) is based on the growth and merging of
nearly empty spheres, motivated by the expectation that voids are spherical to
first order. V2 (Voronoi Voids) is a Python 3 implementation of the ZOBOV
(ZOnes Bordering On Voidness; Neyrinck, 2008) algorithm, which defines voids
based on local minima in the galactic density field. An example of the results
from the VAST void-finding algorithms applied to SDSS DR7 are shown in Figure
1. Any galaxy which falls within one of the void regions defined by a given
algorithm is considered a void galaxy; all others are classified as wall
galaxies.

###  4.1 VoidFinder

VoidFinder begins by identifying and removing all field galaxies, defined as
those galaxies with their third nearest neighbor farther than 5.33$h^{-1}$
Mpc. The remaining (wall) galaxies are placed on a three dimensional grid with
side lengths of 5$h^{-1}$ Mpc. Spheres are grown from the center of each empty
grid cell and expand until they are bounded by four galaxies on their surface.
Each void is defined as a union of these spheres, where the largest sphere
within a void is that void’s maximal sphere. Maximal spheres are not permitted
to overlap each other by more than 10% of their volume, and they have a
minimum radius of 10$h^{-1}$ Mpc. All other spheres of a given void must
overlap only their void’s maximal sphere by at least 50% of their volumes. For
more details, see Section A.1 and Douglass et al. (2023).

###  4.2 Voronoi Voids (V2)

V2 first produces a 3D Voronoi tessellation of the input galaxy catalog. The
Voronoi cells are combined into “zones” using a watershed segmentation. Zones
are then merged into voids if the least-dense pair of adjacent cells between
two zones is below a given maximum “linking density.” Here, we define the
zone-linking density threshold as $0.2\bar{\rho}$, where $\bar{\rho}$ is the
mean density of the SDSS DR7 volume-limited galaxy catalog, similar to the
threshold used in VIDE (Sutter et al., 2015). In addition, we prune our V2
catalog to remove all voids with an effective radius smaller than 10$h^{-1}$
Mpc. For more details, see Section A.2 and Douglass et al. (2023).

Unlike VoidFinder, which assumes that the voids are spherical to first order,
V2 makes no assumptions regarding void shape. In addition, V2 depends on fewer
tunable parameters than VoidFinder, making it more survey-agnostic.

##  5 Comparison of void galaxy properties between algorithms

Table 1: Galaxy property distribution statistics

Algorithm | Environment | Average | Median | Average shift | Median shift | $\log_{10}{(B_{12})}$
---|---|---|---|---|---|---
Stellar mass [$\log(M_{*}/M_{\odot})$]
VoidFinder | Void | 9.786 ±0.002 | 9.910 | $0.205\pm 0.002\text{ (}2.1\pm 0.2$%) | $0.164\text{ (}1.6$%) | $-1708$
Wall | 9.992 ±0.001 | 10.074
V2 | Void | 9.927 ±0.001 | 10.025 | $0.038\pm 0.002\text{ (}0.4\pm 0.2$%) | $0.043\text{ (}0.4$%) | $-203$
Wall | 9.965 ±0.002 | 10.068
Absolute magnitude, $M_{r}$
VoidFinder | Void | -19.619 ±0.004 | -19.837 | $-0.359\pm 0.005\text{ (}1.8\pm 0.5$%) | $-0.263\text{ (}1.3$%) | $-1248$
Wall | -19.978 ±0.002 | -20.100
V2 | Void | -19.859 ±0.002 | -20.013 | $-0.075\pm 0.005\text{ (}0.4\pm 0.5$%) | $-0.088\text{ (}0.4$%) | $-308$
Wall | -19.934 ±0.004 | -20.101
$u-r$
VoidFinder | Void | 1.763 ±0.002 | 1.730 | $0.173\pm 0.002\text{ (}9.0\pm 0.2$%) | $0.259\text{ (}13$%) | $-1592$
Wall | 1.937 ±0.001 | 1.989
V2 | Void | 2.102 ±0.001 | 2.127 | $0.027\pm 0.002\text{ (}1.3\pm 0.2$%) | $0.039\text{ (}1.8$%) | $-28.8$
Wall | 2.130 ±0.002 | 2.166
$g-r$
VoidFinder | Void | 0.5614 ±0.0007 | 0.5729 | $0.0540\pm 0.0008\text{ (}8.77\pm 0.08$%) | $0.0855\text{ (}13.0$%) | $-1704$
Wall | 0.6154 ±0.0004 | 0.6584
V2 | Void | 0.6836 ±0.0005 | 0.7201 | $0.0128\pm 0.0009\text{ (}1.84\pm 0.09$%) | $0.0172\text{ (}2.33$%) | $-57.9$
Wall | 0.6963 ±0.0007 | 0.7373
SFR [$\log(M_{\odot}/\text{yr})$]
VoidFinder | Void | -0.303 ±0.002 | -0.216 | $-0.097\pm 0.003\text{ (}24.2\pm 0.3$%) | $-0.175\text{ (}44.8$%) | $-416$
Wall | -0.400 ±0.001 | -0.391
V2 | Void | -0.385 ±0.001 | -0.355 | $0.008\pm 0.003\text{ (}2.1\pm 0.3$%) | $0.020\text{ (}5.9$%) | $-66.5$
Wall | -0.377 ±0.002 | -0.336
sSFR [$\log(\text{yr}^{-1})$]
VoidFinder | Void | -10.489 ±0.003 | -10.232 | $-0.294\pm 0.003\text{ (}2.7\pm 0.3$%) | $-0.445\text{ (}4.2$%) | $-1471$
Wall | -10.783 ±0.002 | -10.677
V2 | Void | -10.705 ±0.002 | -10.521 | $-0.044\pm 0.003\text{ (}0.4\pm 0.3$%) | $-0.094\text{ (}0.9$%) | $-18.9$
Wall | -10.750 ±0.003 | -10.615

Note. — Galaxy property distribution summary for void and wall environments
according to VoidFinder or V2. The shifts are calculated as $($wall $-$
void$)$, and the percentages are equal to
$(\text{wall}-\text{void})/\text{wall}$, with a negative sign indicating that
void galaxies have larger values than the wall galaxies. The logarithm of the
Bayes factor is included; all Bayes factors are significantly less than
$\log(1/100)=-2$, indicating definitive evidence for the two-parent model
(Appendix A.3).

We investigate the influence of the void environment on inferred galaxy
properties, where the environmental classification (void or wall) is
determined from each of the void catalogs described in Section 4. To study how
the inferred properties of void galaxies depend on the choice of void-finding
algorithm, we plot the distribution of each property for the void and wall
galaxy populations as defined by VoidFinder and V2. We report the mean,
median, and corresponding difference between the means and medians of the void
and wall galaxy distributions for each property and algorithm in Table 1.

Two Bayesian analysis models are used to quantify the statistical significance
of differences between the inferred properties of void and wall galaxies.
Properties with unimodal distributions (luminosity and stellar mass) are
modeled as a mixture of two skew normal distributions; properties with bimodal
distributions (color and (s)SFR) are modeled as a mixture of three skew normal
distributions. We then compute the Poisson likelihoods that the void and wall
galaxy properties are drawn from the same parent distribution and from
different parent populations. The likelihood ratio of the one- and two-parent
models is used to construct a Bayes factor, $B_{12}$, to discriminate between
the models (see Appendix A.3 for details).

For all properties and void-finding algorithms, the Bayes factor gives
decisive evidence in favor of the two-parent model. However, the differences
between void and wall galaxies are substantially less significant for V2 than
VoidFinder. The Bayes factors for all galaxy properties and void-finding
algorithms are listed in the rightmost column of Table 1. In the following
subsections, we describe differences between the stellar mass, color,
luminosity, and (specific) star formation rate of void and wall galaxies using
V2 and VoidFinder for void galaxy classification.

###  5.1 Stellar Mass

We study the distribution of stellar masses of void galaxies compared to the
distribution of galaxies in denser regions, as shown in Figure 2. Void
galaxies are expected to have systematically lower stellar masses than
galaxies in denser regions because galaxy interactions, much more rare in
voids, are expected to enhance star formation. Both the distributions shown in
Figure 2 and the statistics reported in Table 1 show that the stellar mass
distribution for void galaxies is statistically different than the wall galaxy
sample when either void-finding algorithm is used for void classification.
When the void environment is defined with VoidFinder, we see this difference
manifest as a relatively large shift towards lower stellar masses in the
voids. With V2, though, we find that the difference between the void and wall
populations is very slight and appears to be a result of a slight narrowing of
the stellar mass distribution in walls, as the void galaxy distribution
appears to match the distribution of all galaxies.

Figure 2: Stellar mass distribution of galaxies, separated into void and wall
environments according to the VoidFinder (left) and V2 (right) void catalogs.
The distribution of all galaxies is shaded in gray. There is a shift towards
lower stellar mass for void galaxies relative to wall galaxies with both void
classifications, although V2’s is very slight.

###  5.2 Luminosity

Figure 3: Luminosity distribution of galaxies separated according to void and
wall environments by VoidFinder (left) and V2 (right). The distribution for
all galaxies is shaded in gray. There is a shift towards fainter luminosities
in the distribution of void galaxies relative to wall galaxies in VoidFinder’s
classification, while it appears that the wall galaxies are shifted towards
slightly brighter luminosities with V2’s classification.

A galaxy’s luminosity is closely linked to its stellar mass, so we expect void
galaxies to be systematically fainter than galaxies in denser regions. As seen
on the left in Figure 3, the classification by VoidFinder results in a shift
towards fainter luminosities for void galaxies than galaxies in denser
environments. This matches the trend observed in the stellar mass, shown in
Figure 2. Using the V2 void catalog, though, we see very little difference
between the void and wall galaxy luminosity distributions, as seen on the
right in Figure 3; the void galaxies as defined by V2 have nearly the same
luminosities as its wall galaxies, with possibly a slight thinning of the wall
population. The Bayes factors (listed in Table 1) corroborate these
differences between the void and wall populations for both classifications.

###  5.3 (specific) Star formation rate

Figure 4: Star formation rate (SFR, top) and specific star formation rate
(sSFR, bottom) distributions of galaxies for void and wall environments as
classified by VoidFinder (left) and V2 (right). The distribution of all
galaxies is shown in gray for comparison. When classified using VoidFinder’s
void catalog, void galaxies have systematically higher (specific) star
formation rates relative to wall galaxies; we observe no discernible
difference between void and wall galaxies when the V2 void catalog is used.

Void galaxies are expected to be undergoing more star formation today than
galaxies in denser regions, due both to potential retarded star formation and
their continued presence of cold gas reservoirs (Cen, 2011). Therefore, void
galaxies are expected to have higher star formation rates (SFR) and specific
star formation rates (sSFR; the star formation rate per unit mass) than
galaxies in denser regions. As shown in the left column of Figure 4, we find
that void galaxies have higher (s)SFR values relative to galaxies in denser
environments when using the VoidFinder void catalog for environmental
classification. As seen on the right in Figure 4, the distributions for V2
show no such trend in (s)SFR — according to the V2 void catalog, both void and
wall galaxies have the same (s)SFR.

The Bayes factors (listed in Table 1) support the significant difference
between the void and wall populations of VoidFinder for the (s)SFR
distributions. However, the Bayes factors for the (s)SFR distributions for
void and wall galaxies as classified by V2 also indicate that the two
populations come from different parent distributions. From the distributions
shown on the right in Figure 4, this is possibly due to a very slight shift in
the V2 wall galaxy population towards lower (s)SFRs, not a shift in the void
galaxy population towards higher (s)SFRs as seen with VoidFinder.

###  5.4 Color

Figure 5: Color distribution, $u-r$ (top row) and $g-r$ (bottom row),
separated by their environment (void, wall) using VoidFinder (left column) and
V2 (right column). The distribution for all galaxies is given in gray. There
is an apparent shift towards bluer color for void galaxies relative to wall
galaxies with VoidFinder’s classification but not V2’s.

For galaxies in SDSS DR7, we explore the distribution of the color of void
galaxies relative to galaxies in denser regions as classified by both
VoidFinder and V2, as shown in Figure 5. A galaxy’s color is determined by its
stellar populations: bluer colors are a result of hot O and B stars,
indicating that star formation has recently occurred, and red colors are
indicative of an old stellar population. With star formation expected to shift
to later times in void regions, void galaxies are expected to be
systematically bluer than galaxies in denser regions. This is exactly the
shift that we see in the VoidFinder void galaxy population, shown on the left
in Figure 5 and reflected in the Bayes factor reported in Table 1 comparing
the VoidFinder galaxy distributions.

The Bayes factors for the V2 color distributions are two orders of magnitude
smaller than those of VoidFinder. Combined with the distributions shown on the
right in Figure 5, we conclude that the differences quantified by the Bayes
factors in the V2 distributions are actually a result of a slight increase in
redder galaxies in the wall population, as evidenced by the slight increase in
amplitude of the right-hand peak and slight decrease in the left-hand peak of
the V2 wall distributions compared to the void and overall distributions.

##  6 Discussion

In Section 5, we compare the distributions of a range of galaxy properties for
galaxies in voids and galaxies in denser regions using two different
definitions of the void environment. With VoidFinder, we find that void
galaxies are systematically fainter, bluer, less massive, and have
considerably higher star formation rates than galaxies in denser regions. With
V2, though, we see very little difference between the properties of void
galaxies and the properties of galaxies in denser regions.

###  6.1 Impact of the void-finding algorithm on galaxy classification

The $\Lambda$CDM cosmology predicts void galaxies to be retarded in their star
formation relative to galaxies in denser regions since they experience fewer
interactions (which often initiate star formation episodes) compared to
galaxies in denser regions. Additionally, they are able to retain their gas as
they evolve in these extremely underdense environments (Cen, 2011), allowing
them to continue forming stars up to the present epoch. The void galaxy
properties with VoidFinder used to define void regions are therefore
consistent with the simulated physical properties of void galaxies (Benson et
al., 2003; Cen, 2011; Habouzit et al., 2020; Peper & Roukema, 2021) and with
previous observations: void galaxies are fainter (e.g, Croton et al., 2005;
Hoyle et al., 2005; Sorrentino et al., 2006; Kreckel et al., 2012; Moorman et
al., 2015), less massive (e.g., Jung et al., 2014; Alpaslan et al., 2016;
Beygu et al., 2017; Martizzi et al., 2020), bluer (e.g., Grogin & Geller,
1999; Rojas et al., 2004; Croton et al., 2005; Patiri et al., 2006; Sorrentino
et al., 2006; von Benda-Beckmann & Müller, 2008; Hoyle et al., 2012; Kreckel
et al., 2012), and have higher (s)SFR (e.g., Rojas et al., 2005; von Benda-
Beckmann & Müller, 2008; Cybulski et al., 2014; Liu et al., 2015; Beygu et
al., 2016; Moorman et al., 2016).

However, we do not see these expected traits for galaxies in voids defined by
V2. To better understand the source of this discrepancy, we investigate the V2
void galaxies in more detail. We find that 23.4% of VoidFinder wall galaxies
are classified as void galaxies by V2, and 40.5% of V2 void galaxies are
classified as wall galaxies by VoidFinder. As described in Section 4,
VoidFinder voids are found by expanding spheres until they reach a relatively
bright galaxy that is not isolated, while V2 defines voids by connecting
Voronoi cells until they reach local density maxima. Since these density
maxima are likely to be close to the centers of the walls/filaments rather
than their edges (where the brighter non-isolated galaxies that VoidFinder
stops at are located), V2 likely extends further into what VoidFinder
considers a wall. This leads to V2 classifying many galaxies as void that are,
instead, classified as wall galaxies by VoidFinder.

Considering how the algorithms define the edges of voids and how VoidFinder’s
results are consistent with expectations, we conclude that the void
environment of V2 is contaminated by wall galaxies, similar to that shown by
Florez et al. (2021) and Veyrat et al. (2023). V2’s void environment is
therefore expected to behave more like the general distribution of galaxies.
Our observations agree with the simulation results of Veyrat et al. (2023):
VoidFinder provides a better delineation between the dynamically-distinct void
regions and the surrounding wall structures, leading to more reasonable
physical conclusions regarding the properties of void galaxies. Unless
additional restrictions are implemented, this contamination results in the
void-wall classification from V2, and thus ZOBOV-like algorithms, to not be
suitable for void environment studies of galaxy formation and evolution.

###  6.2 Galaxy properties as a function of depth within V2 voids

Figure 6: The relationship between the distribution of sSFR of void galaxies,
as classified by V2, and the galaxies’ depths within the voids in which they
reside (whether they live close to the edge of the void or deep inside it).
The deepest galaxies within each void have a normalized distance of 1, while a
distance closer to 0 indicates galaxies that are close to the edge. Galaxies
further from the edges of the V2 voids have sSFR distributions which more
closely match that of VoidFinder void galaxies; as galaxies closer to the V2
void boundaries are included in the sample, the sSFR distribution shifts to
more closely resemble the SDSS DR7 galaxy distribution.

If the void samples of V2 are contaminated by wall galaxies located at the
edges of its voids, the expected void galaxy properties can be recovered by
focusing on the galaxies that are closer to the centers of its voids. Studying
the galaxy properties as a function of void depth, we can determine the void
depth within which V2 void galaxies have similar trends to VoidFinder void
galaxies. This is the observational complement to one of the mitigation
techniques discussed in Veyrat et al. (2023), who tested this idea on the
dynamical classification of simulated dark matter particles.

We calculate a galaxy’s distance to the void’s edge to determine how deep a
galaxy is within a V2 void. This is equivalent to the distance from the
position of the galaxy to the edge of the void, defined by the Voronoi-cell
surfaces that make up the boundary between a cell inside the void and a cell
outside of it. We normalize the depth of each galaxy by the largest depth
calculated for that void, so that we can compare voids of different sizes.
Therefore, the deepest galaxy within each void has a depth of 1, and a galaxy
close to the edge will have a depth close to 0.

To test this hypothesis, we fit subsets of the V2 void galaxies with a skew
normal (or sum of two skew normals if the distribution is bimodal) as a
function of void depth for each galaxy property studied in Section 5. For a
given void depth, the corresponding subsample of void galaxies consists of all
galaxies with depths greater than the given depth threshold. This allows us to
study the behavior of the V2 void galaxy population if the voids were limited
to some inner fraction of the current total void size. The results of this
analysis for the sSFR is shown in Figure 6. For comparison, we also include
the ratio of the sSFR peak amplitudes for VoidFinder void galaxies (in red)
and the ratio of the sSFR peak amplitudes for all galaxies in the overall NSA
catalog (in black).

For unimodal galaxy properties (stellar mass and luminosity), we find a steady
shift from smaller masses / fainter luminosities to larger masses / brighter
luminosities as the V2 void galaxies include more of those closer to the void
boundaries. For properties with bimodal distributions (color and (s)SFR), the
fraction of blue galaxies / galaxies with higher (s)SFR decreases while the
fraction of red galaxies / galaxies with lower (s)SFR increases as we include
more V2 void galaxies closer to their void boundaries. Figure 6 shows the
ratio of the peak of the mode at higher sSFR to the mode at lower sSFR as a
function of minimum depth included in the subsample.

For all properties, we find that the V2 void galaxy distribution matches that
of VoidFinder void galaxies when only V2 void galaxies in the inner $\sim$20%
of the volume of the V2 voids are included in the void sample. Therefore, the
galaxies sitting close to the boundaries of V2 voids are more likely to be
wall galaxies and not actual void galaxies. This matches the results of
Habouzit et al. (2020), who apply VIDE to a cosmological hydrodynamical
simulation and reproduce the expected void galaxy property results when
including only those galaxies closest to the centers of the VIDE voids. This
also matches the simulation results of Veyrat et al. (2023), who compare the
distribution of dark matter particles with a given crossing number as a
function of void depth for both VoidFinder and V2. Nadathur & Hotchkiss (2014)
also discusses the tendency for the watershed algorithm to include high-
density walls and filaments within void regions.

Therefore, our analysis of galaxy properties as a function of depth within V2
voids is consistent with our hypothesis and previous studies that the void
sample of V2 is contaminated by wall galaxies, a result of the extension of V2
voids into the walls. This is visible on the right in Figure 1, where
clustered void galaxies (in red) are located at the void edges. This
phenomenon is also visible in our mock 2D galaxy catalog (see the middle and
right panels of Figure 9), where multiple cells outside the artificial voids
are joined to the void region. This “leaking” of the void regions into the
walls is a consequence of how the watershed algorithm within V2 (and other
ZOBOV-like void-finders) defines void boundaries based on density maxima. To
use V2 voids to study galaxy properties, one must therefore implement a depth
limit to exclude wall galaxies from the void galaxy population.

##  7 Conclusions

We examine the impact of the choice of void-finding algorithm used to classify
the large-scale structure of galaxies by comparing properties of galaxies
determined to be void galaxies using the VoidFinder and V2 algorithms. Using
the VoidFinder void catalog for SDSS DR7, we find that the distributions of
void galaxy properties are consistent with the expected physical properties of
galaxies existing in voids: smaller masses, fainter luminosities, bluer
colors, and higher (specific) star formation rates. When using V2, however, we
find very little difference between the properties of void galaxies and the
properties of galaxies in denser regions.

Upon further investigation, we find that V2 void regions include a significant
fraction of the walls and filaments. This is a result of the watershed
algorithm extending void regions to the galaxy density maxima, which are
located in the centers of the denser regions in the galaxy distribution. While
VoidFinder is better suited for studying void galaxy properties because it is
not subject to this phenomenon, it is possible to recover the void galaxy
population with V2 voids if voids are limited to the inner $\sim$20% of their
volume.

##  8 Acknowledgements

The authors thank Stephen W. O’Neill, Jr. for his help and expertise with
VAST, as well as Michael Vogeley, Regina Demina, and Tom Ferbel for their
feedback on this manuscript. D. Veyrat and S. BenZvi acknowledge support from
the U.S. Department of Energy Office of High Energy Physics under the Award
Number DE-SC0008475.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan
Foundation, the Participating Institutions, the National Science Foundation,
the U.S. Department of Energy, the National Aeronautics and Space
Administration, the Japanese Monbukagakusho, the Max Planck Society, and the
Higher Education Funding Council for England. The SDSS web site is
http://www.sdss.org/.

The SDSS is managed by the Astrophysical Consortium for the Participating
Institutions. The Participating Institutions are the American Museum of
Natural History, Astrophysical Institute Potsdam, University of Basel,
University of Cambridge, Case Western Reserve University, University of
Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the
Japan Participation Group, Johns Hopkins University, the Joint Institute for
Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and
Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences
(LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for
Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico
State University, Ohio State University, University of Pittsburgh, University
of Portsmouth, Princeton University, the United States Naval Observatory, and
the University of Washington.

##  Appendix A VAST examples

Figure 7: A sample random 2D galaxy catalog with circles of radius $10\bar{r}$
removed in each quadrant to simulate simple voids, where $\bar{r}$ is the mean
inter-particle separation.

To visualize how the two algorithms included in VAST identify voids, we apply
each algorithm to a two-dimensional (2D) mock catalog of randomly distributed
galaxies. The artificial voids are formed by emptying four circular holes in
the random galaxy distribution, as shown in Figure 7.

###  A.1 VoidFinder example

Figure 8: Steps in the VoidFinder algorithm applied to a 2D mock galaxy
catalog. _Left:_ Non-field (wall) galaxies placed on a grid. Any grid cell
containing at least one galaxy is shaded red; all empty grid cells are blue.
_Center:_ Circles grown from the center of each empty grid cell until bounded
by at least three galaxies. _Right:_ Voids formed through the union of the
grown spheres.

After removing isolated (field) galaxies, the remaining (wall) galaxies are
placed on a 5$h^{-1}$ Mpc grid. An example of this grid is shown on the left
in Figure 8 for the 2D mock galaxies, where empty cells are colored blue.

Next, a sphere (or hole) is expanded from each empty cell in the grid. Each
sphere, centered on the cell, is created with radius equal to the distance to
the nearest wall galaxy. The sphere’s center then shifts away from that galaxy
as its radius increases until its surface encounters a second galaxy. At this
point, the sphere expands from the midpoint between these two galaxies until a
third galaxy is encountered. Finally, the sphere grows perpendicularly away
from the plane formed by these three galaxies until it reaches a fourth
galaxy. The center panel of Figure 8 shows the result of this step for the 2D
catalog of mock galaxies.

To extract unique voids from this list of holes, the holes are sorted by
radius, and the largest hole is identified as a maximal sphere — the largest
sphere that can fit in a given void region. Then, for each subsequent hole, if
it does not overlap any already-identified maximal spheres by more than 10% of
its volume, it is also identified as a maximal sphere. This process continues
for all holes with radii larger than $10\bar{r}$, where $\bar{r}$ is the mean
particle separation in the simulated random sample.

Finally, the volume of the void is refined by merging each maximal sphere with
the remaining non-maximal spheres. When a hole overlaps only one maximal
sphere by more than 50% of its volume, it is merged into that void. Holes that
overlap multiple maximal spheres by more than 50% are discarded. Each void is
therefore composed of one maximal sphere and some number of smaller holes. The
right panel in Figure 8 shows the final 2D catalog of voids found using
VoidFinder.

###  A.2 V2 example

Figure 9: Steps in the V2 algorithm applied to a 2D mock galaxy catalog.
_Left:_ The Voronoi tessellation of the catalog. _Center:_ Zones formed using
a watershed method. _Right:_ Final voids identified after merging zones and
pruning voids.

V2 can be divided into four steps:

  1. 1.

Tessellation: A Voronoi tessellation of the input catalog is produced. The
tessellation of the 2D mock galaxies is shown on the left in Figure 9. The
inverse of each cell’s volume is an estimate of that galaxy’s local density.

  2. 2.

Zone formation: Voronoi cells are combined into groups or “zones” using
watershed segmentation. Each cell is put into the same zone as its least-dense
neighbor, and cells less dense than any neighboring cell (local density
minima) are identified as their zone’s central cell. The center panel of
Figure 9 shows the zones formed using the 2D mock sample.

  3. 3.

Zone merging: Zones are merged to form voids by first identifying the least-
dense pair of adjacent cells between two zones. For a given maximum “linking
density,” the set of all zones is partitioned into a subset of voids linked
together by zones of equal or lower density. The algorithm then loops over all
linking densities to produce a list of voids.

  4. 4.

Void pruning: Zone linking produces a hierarchy of voids ranging from
individual zones to a single void encompassing the entire survey volume. A
pruning step similar to that implemented in VIDE (Sutter et al., 2015) is
therefore introduced to identify the final list of voids. The final result for
the 2D sample with this pruning is shown on the right in Figure 9.

###  A.3 Calculating the Bayes factor

Due to our large sample size, non-Bayesian significance tests such as a two-
sample Kolmogorov–Smirnov (K-S) test yield extremely low p-values that
inaccurately quantify the extent of similarity of our distributions. Since we
know that our data is distributed like one or more Gaussian or skew normal
distributions (a mixture model), we can bin and fit the data to the model
assuming that the samples come from either the same parent distribution or two
different parent distributions. Then, we can compare these two models using
the Bayes Factor.

We have two data sets with binned counts $\mathbf{m}={m_{i}}$ and
$\mathbf{n}={n_{i}}$. For the single-parent model, we assume that both data
sets are created by the same skew normal mixture model, such that the mean
count in bin $i$ for data set 1 is

$\lambda_{1,i}=\left[\alpha\
\mathrm{SN}(\mu_{\alpha},\sigma_{\alpha},\xi_{\alpha})+\beta\
\mathrm{SN}(\mu_{\beta},\sigma_{\beta},\xi_{\beta})\right]_{i}.$ (A1)

To account for the different sizes in the data sets, we can add a scale factor
$s$ to the fit of the second data set:

$\lambda_{2,i}=\left[s\left(\alpha\
\mathrm{SN}(\mu_{\alpha},\sigma_{\alpha},\xi_{\alpha})+\beta\
\mathrm{SN}(\mu_{\beta},\sigma_{\beta},\xi_{\beta})\right)\right]_{i}=s\lambda_{1,i}.$
(A2)

Thus the model includes 9 parameters
$\theta=\left\\{s,\alpha,\mu_{\alpha},\sigma_{\alpha},\xi_{\alpha},\beta,\mu_{\beta},\sigma_{\beta},\xi_{\beta}\right\\}$.
We can then maximize the joint Poisson likelihood

$p(\mathbf{m},\mathbf{n}|\lambda_{1},\lambda_{2},\mathcal{M}_{1})=\prod_{i=1}^{N}\frac{\lambda_{1,i}^{m_{i}}e^{-\lambda_{1,i}}}{m_{i}!}\frac{\lambda_{2,i}^{n_{i}}e^{-\lambda_{2,i}}}{n_{i}!},$
(A3)

or minimize the negative log likelihood

$-\ln{\mathcal{L}}=-\sum_{i=1}^{N}\left(m_{i}\lambda_{1,i}-\lambda_{1,i}-\ln{m_{i}!}+n_{i}\lambda_{2,i}-\lambda_{2,i}-\ln{n_{i}!}\right)$
(A4)

to obtain the best-fit parameters.

If on the other hand we assume the binned counts $\mathbf{m}={m_{i}}$ and
$\mathbf{n}={n_{i}}$ are created by separate skew normal mixture models with
mean count

$\lambda_{1,i}=\left[\alpha\
\mathrm{SN}(\mu_{\alpha},\sigma_{\alpha},\xi_{\alpha})+\beta\
\mathrm{SN}(\mu_{\beta},\sigma_{\beta},\xi_{\beta})\right]_{i}$ (A5)

in set 1 and mean count

$\lambda_{2,i}=\left[\gamma\
\mathrm{SN}(\mu_{\gamma},\sigma_{\gamma},\xi_{\gamma})+\delta\
\mathrm{SN}(\mu_{\delta},\sigma_{\delta},\xi_{\delta})\right]_{i}$ (A6)

in set 2, then the model includes 16 free parameters

$\theta=\left\\{\alpha,\mu_{\alpha},\sigma_{\alpha},\xi_{\alpha},\beta,\mu_{\beta},\sigma_{\beta},\xi_{\beta},\gamma,\mu_{\gamma},\sigma_{\gamma},\xi_{\gamma},\delta,\mu_{\delta},\sigma_{\delta},\xi_{\delta}\right\\}.$

We can still use the likelihood in Eqn. A3 with mean counts from Eqns. A5 and
A6.

To evaluate the support of the Bayes factor for the one or two-parent
hypothesis, we use the thresholds suggested in Jeffreys (1961) and Kass &
Raftery (1995).

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543, doi: 10.1088/0067-0049/182/2/543
  * Alpaslan et al. (2016) Alpaslan, M., Grootes, M., Marcum, P. M., et al. 2016, MNRAS, 457, 2287, doi: 10.1093/mnras/stw134
  * Benson et al. (2003) Benson, A. J., Hoyle, F., Torres, F., & Vogeley, M. S. 2003, MNRAS, 340, 160, doi: 10.1046/j.1365-8711.2003.06281.x
  * Beygu et al. (2016) Beygu, B., Kreckel, K., van der Hulst, J. M., et al. 2016, MNRAS, 458, 394, doi: 10.1093/mnras/stw280
  * Beygu et al. (2017) Beygu, B., Peletier, R. F., van der Hulst, J. M., et al. 2017, MNRAS, 464, 666, doi: 10.1093/mnras/stw2362
  * Blanton et al. (2011) Blanton, M. R., Kazin, E., Muna, D., Weaver, B. A., & Price-Whelan, A. 2011, AJ, 142, 31, doi: 10.1088/0004-6256/142/1/31
  * Blanton et al. (2003) Blanton, M. R., Lin, H., Lupton, R. H., et al. 2003, AJ, 125, 2276, doi: 10.1086/344761
  * Bond et al. (1996) Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603, doi: 10.1038/380603a0
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Cen (2011) Cen, R. 2011, ApJ, 741, 99, doi: 10.1088/0004-637X/741/2/99
  * Colberg et al. (2008) Colberg, J. M., Pearce, F., Foster, C., et al. 2008, MNRAS, 387, 933, doi: 10.1111/j.1365-2966.2008.13307.x
  * Constantin et al. (2008) Constantin, A., Hoyle, F., & Vogeley, M. S. 2008, ApJ, 673, 715, doi: 10.1086/524310
  * Croton et al. (2005) Croton, D. J., Farrar, G. R., Norberg, P., et al. 2005, MNRAS, 356, 1155, doi: 10.1111/j.1365-2966.2004.08546.x
  * Cybulski et al. (2014) Cybulski, R., Yun, M. S., Fazio, G. G., & Gutermuth, R. A. 2014, MNRAS, 439, 3564, doi: 10.1093/mnras/stu200
  * de Lapparent et al. (1986) de Lapparent, V., Geller, M. J., & Huchra, J. P. 1986, ApJ, 302, L1, doi: 10.1086/184625
  * Douglass et al. (2022) Douglass, K., Veyrat, D., O’Neill, S., et al. 2022, JOSS, 7, 4033, doi: 10.21105/joss.04033
  * Douglass et al. (2023) Douglass, K. A., Veyrat, D., & BenZvi, S. 2023, ApJS, 265, 7, doi: 10.3847/1538-4365/acabcf
  * Douglass & Vogeley (2017) Douglass, K. A., & Vogeley, M. S. 2017, ApJ, 837, 42, doi: 10.3847/1538-4357/aa5e53
  * El-Ad & Piran (1997) El-Ad, H., & Piran, T. 1997, ApJ, 491, 421, doi: 10.1086/304973
  * Florez et al. (2021) Florez, J., Berlind, A. A., Kannappan, S. J., et al. 2021, ApJ, 906, 97, doi: 10.3847/1538-4357/abca9f
  * Grogin & Geller (1999) Grogin, N. A., & Geller, M. J. 1999, AJ, 118, 2561, doi: 10.1086/301126
  * Grogin & Geller (2000) —. 2000, AJ, 119, 32, doi: 10.1086/301179
  * Habouzit et al. (2020) Habouzit, M., Pisani, A., Goulding, A., et al. 2020, MNRAS, 493, 899, doi: 10.1093/mnras/staa219
  * Hamaus et al. (2020) Hamaus, N., Pisani, A., Choi, J.-A., et al. 2020, J. Cosmology Astropart. Phys, 2020, 023, doi: 10.1088/1475-7516/2020/12/023
  * Hamaus et al. (2016) Hamaus, N., Pisani, A., Sutter, P. M., et al. 2016, Phys. Rev. Lett., 117, 091302, doi: 10.1103/PhysRevLett.117.091302
  * Hoyle et al. (2005) Hoyle, F., Rojas, R. R., Vogeley, M. S., & Brinkmann, J. 2005, ApJ, 620, 618, doi: 10.1086/427176
  * Hoyle et al. (2012) Hoyle, F., Vogeley, M. S., & Pan, D. 2012, MNRAS, 426, 3041, doi: 10.1111/j.1365-2966.2012.21943.x
  * Jeffreys (1961) Jeffreys, H. 1961, Theory of Probability, 3rd edn. (Oxford, UK: Oxford University Press)
  * Jung et al. (2014) Jung, I., Lee, J., & Yi, S. K. 2014, ApJ, 794, 74, doi: 10.1088/0004-637X/794/1/74
  * Kass & Raftery (1995) Kass, R. E., & Raftery, A. E. 1995, Journal of the American Statistical Association, 90, 773, doi: 10.1080/01621459.1995.10476572
  * Kreckel et al. (2012) Kreckel, K., Platen, E., Aragón-Calvo, M. A., et al. 2012, AJ, 144, 16, doi: 10.1088/0004-6256/144/1/16
  * Lavaux & Wandelt (2012) Lavaux, G., & Wandelt, B. D. 2012, ApJ, 754, 109, doi: 10.1088/0004-637X/754/2/109
  * Liu et al. (2015) Liu, C.-X., Pan, D. C., Hao, L., et al. 2015, ApJ, 810, 165, doi: 10.1088/0004-637X/810/2/165
  * Lupton et al. (2001) Lupton, R., Gunn, J. E., Ivezić, Z., Knapp, G. R., & Kent, S. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 238, Astronomical Data Analysis Software and Systems X, ed. F. R. Harnden, Jr., F. A. Primini, & H. E. Payne, 269
  * Mao et al. (2017) Mao, Q., Berlind, A. A., Scherrer, R. J., et al. 2017, ApJ, 835, 160, doi: 10.3847/1538-4357/835/2/160
  * Martizzi et al. (2020) Martizzi, D., Vogelsberger, M., Torrey, P., et al. 2020, MNRAS, 491, 5747, doi: 10.1093/mnras/stz3418
  * Moorman et al. (2016) Moorman, C. M., Moreno, J., White, A., et al. 2016, ApJ, 831, 118, doi: 10.3847/0004-637X/831/2/118
  * Moorman et al. (2015) Moorman, C. M., Vogeley, M. S., Hoyle, F., et al. 2015, ApJ, 810, 108, doi: 10.1088/0004-637X/810/2/108
  * Muldrew et al. (2012) Muldrew, S. I., Croton, D. J., Skibba, R. A., et al. 2012, MNRAS, 419, 2670, doi: 10.1111/j.1365-2966.2011.19922.x
  * Nadathur et al. (2019) Nadathur, S., Carter, P. M., Percival, W. J., Winther, H. A., & Bautista, J. E. 2019, Phys. Rev. D, 100, 023504, doi: 10.1103/PhysRevD.100.023504
  * Nadathur & Hotchkiss (2014) Nadathur, S., & Hotchkiss, S. 2014, MNRAS, 440, 1248, doi: 10.1093/mnras/stu349
  * Neyrinck (2008) Neyrinck, M. C. 2008, MNRAS, 386, 2101, doi: 10.1111/j.1365-2966.2008.13180.x
  * Park et al. (2007) Park, C., Choi, Y.-Y., Vogeley, M. S., et al. 2007, ApJ, 658, 898, doi: 10.1086/511059
  * Patiri et al. (2006) Patiri, S. G., Prada, F., Holtzman, J., Klypin, A., & Betancort-Rijo, J. 2006, MNRAS, 372, 1710, doi: 10.1111/j.1365-2966.2006.10975.x
  * Peebles (2001) Peebles, P. J. E. 2001, ApJ, 557, 495, doi: 10.1086/322254
  * Peper & Roukema (2021) Peper, M., & Roukema, B. F. 2021, MNRAS, 505, 1223, doi: 10.1093/mnras/stab1342
  * Pisani et al. (2015) Pisani, A., Sutter, P. M., Hamaus, N., et al. 2015, Phys. Rev. D, 92, 083531, doi: 10.1103/PhysRevD.92.083531
  * Rojas et al. (2004) Rojas, R. R., Vogeley, M. S., Hoyle, F., & Brinkmann, J. 2004, ApJ, 617, 50, doi: 10.1086/425225
  * Rojas et al. (2005) —. 2005, ApJ, 624, 571, doi: 10.1086/428476
  * Sahlén (2019) Sahlén, M. 2019, Phys. Rev. D, 99, 063525, doi: 10.1103/PhysRevD.99.063525
  * Sorrentino et al. (2006) Sorrentino, G., Antonuccio-Delogu, V., & Rifatto, A. 2006, A&A, 460, 673, doi: 10.1051/0004-6361:20065789
  * Strauss et al. (2002) Strauss, M. A., Weinberg, D. H., Lupton, R. H., et al. 2002, AJ, 124, 1810, doi: 10.1086/342343
  * Sutter et al. (2012a) Sutter, P. M., Lavaux, G., Wandelt, B. D., & Weinberg, D. H. 2012a, ApJ, 761, 44, doi: 10.1088/0004-637X/761/1/44
  * Sutter et al. (2012b) —. 2012b, ApJ, 761, 187, doi: 10.1088/0004-637X/761/2/187
  * Sutter et al. (2014) Sutter, P. M., Pisani, A., Wandelt, B. D., & Weinberg, D. H. 2014, MNRAS, 443, 2983, doi: 10.1093/mnras/stu1392
  * Sutter et al. (2015) Sutter, P. M., Lavaux, G., Hamaus, N., et al. 2015, Astronomy and Computing, 9, 1, doi: 10.1016/j.ascom.2014.10.002
  * Veyrat et al. (2023) Veyrat, D., Douglass, K. A., & BenZvi, S. 2023, ApJ, 958, 59, doi: 10.3847/1538-4357/acf4f5
  * von Benda-Beckmann & Müller (2008) von Benda-Beckmann, A. M., & Müller, V. 2008, MNRAS, 384, 1189, doi: 10.1111/j.1365-2966.2007.12789.x
  * Zhao et al. (2020) Zhao, C., Chuang, C.-H., Kitaura, F.-S., et al. 2020, MNRAS, 491, 4554, doi: 10.1093/mnras/stz3339
  * Zhao et al. (2022) Zhao, C., Variu, A., He, M., et al. 2022, MNRAS, 511, 5492, doi: 10.1093/mnras/stac390


11footnotetext: Corresponding author. E-mail<EMAIL_ADDRESS>Fax
(403)-282-5150. Phone (403)-220-3956. Partially supported by the NSERC
Research Grant RGPIN-2020-03934.

# Exponential stability of systems of vector delay differential equations with
applications to second order equations

Leonid Berezansky
Dept. of Math, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel
and Elena Braverman1
Dept. of Math & Stats, University of Calgary, 2500 University Dr. NW, Calgary,
AB, Canada T2N1N4

###### Abstract.

Various results and techniques, such as Bohl-Perron theorem, a priori solution
estimates, M-matrices and the matrix measure, are applied to obtain new
explicit exponential stability conditions for the system of vector functional
differential equations

$\dot{x_{i}}(t)=A_{i}(t)x_{i}(h_{i}(t))+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(t)x_{j}(h_{ij}^{k}(t))+\sum_{j=1}^{n}\int\limits_{g_{ij}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds,~{}i=1,\dots,n.$

Here $x_{i}$ are unknown vector functions, $A_{i},B_{ij}^{k},K_{ij}$ are
matrix functions, $h_{i},h_{ij}^{k},g_{ij}$ are delayed arguments. Using these
results, we deduce explicit exponential stability tests for second order
vector delay differential equations.

Keywords: Exponential stability; differential systems with matrix coefficients
and a distributed delay; second order vector delay differential equations;
Bohl-Perron theorem; matrix measure; M-matrices.

AMS(MOS) subject classification: 34K20, 34K06, 34K25.

##  1\. Introduction

Exponential stability of vector, scalar and higher order delay differential
equations has attracted a lot of attention, see [2, 3, 4, 8, 11, 12, 13, 14,
15, 16] and the bibliography therein. Considered models included neutral
equations and systems with both concentrated and distributed delays. However,
for two classes of linear functional differential equations there are only few
stability results. These types are systems of several linear delay
differential equations of the first order with matrix coefficients and
systems, where at least one equation is of the second or higher order.
Scientific interest to these models is not purely theoretical. There are many
interesting real-world applications of such systems. For example, the ordinary
differential vector equation of the second order

(1.1) $\ddot{x}(t)+A(t)\dot{x}(t)+B(t)x(t)=0$

arises in control theory [9, 17]. Explicit asymptotic stability conditions for
vector ordinary differential equations of the second order were obtained in
[9, 10, 17, 19, 22]. Delay differential vector equations of the second order
are natural generalizations of (1.1) that also can be applied in control
theory. However, presently this application is not practical due to lack of
qualitative results for this class of equations. We can mention only the paper
[20], where for the nonlinear vector equation

(1.2) $\ddot{x}(t)+F(x(t),\dot{x}(t))\dot{x}(t)+H(x(t-\tau))=L(t)$

asymptotic stability of the zero solution and the boundedness of all solutions
of equation (1.2) were investigated. One of the goals of this paper is to fill
this gap, by studying exponential stability of linear and nonlinear vector
delay differential equations of the second order.

The paper is organized as follows. Section 2 contains relevant definitions,
notations and auxiliary statements. In Section 3, we obtain the main result of
the paper on exponential stability for linear systems of functional
differential equations with matrix coefficients. We also give several
corollaries, one of them will further be applied to second order vector
equations. Section 4 deals with exponential stability for a linear delay
differential equation of the second order with matrix coefficients. A
stability test obtained here is one of the first stability results for this
class of equations. As a corollary, we deduce new stability conditions for
non-delay vector equation (1.1). In Section 5, we compare stability results of
the present paper with known ones. Section 5 also includes some open problems
and topics for future research.

##  2\. Preliminaries

The main object of the present paper is

(2.1)
$\dot{x_{i}}(t)=A_{i}(t)x_{i}(h_{i}(t))+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(t)x_{j}(h_{ij}^{k}(t))+\displaystyle\sum_{j=1}^{n}\int\limits_{g_{ij}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds,~{}i=1,\dots,n,~{}t\geq
t_{0}\geq 0,$

where $x_{i}:{\mathbb{R}}\to{\mathbb{R}}^{d}$, $i=1,\dots,n$ are unknown
vector-functions, $A_{i},B_{ij}^{k},K_{ij}$ are given $d\times d$ matrices
with Lebesgue measurable essentially bounded on the semi-axis $[0,\infty)$
entries, and the delays $h_{i},h_{ij}^{k},g_{ij}$ are Lebesgue measurable
functions.

Let us denote by $\|\cdot\|$ an arbitrary vector norm, the same notation will
be used for the induced matrix norm, $I$ is the identity matrix. The absolute
value of the matrix is understood componentwise: if
$B=\left(b_{ij}\right)_{i,j=1}^{d}$ then
$|B|:=\left(|b_{ij}|\right)_{i,j=1}^{d}$, we use the same notation for
vectors. Matrix or vector inequalities $A\leq B$ or $A<B$ are also understood
componentwise: each entry satisfies the inequality. All the relations
(equations and inequalities) with Lebesgue measurable functions are understood
almost everywhere.

To get main results, we will utilize the matrix measure (the logarithmic norm)
defined in [18, 21] and [14, Table 3.1, p. 286] as

$\mu(B):=\lim_{\vartheta\rightarrow 0^{+}}\frac{\|I+\vartheta
B\|-1}{\vartheta}\,.$

In particular, for the maximum norm
${\|x\|}_{\infty}:=\max_{i=1,\dots,d}|x_{i}|$ we have for
$B=(b_{ij})_{i,j=1}^{d}$

${\|B\|}_{\infty}:=\max_{i=1,\dots,d}\sum_{j=1}^{d}|b_{ij}|,~{}~{}\mu_{\infty}(B)=\max_{i=1\dots,d}\left\\{b_{ii}+\sum_{j=1,j\neq
i}^{d}|a_{ij}|\right\\}.$

The following classical definition for an $M$-matrix will be used.

######  Definition 2.1.

[6] A matrix $B=(b_{ij})_{i,j=1}^{n}$ is called a (non-singular) $M$-matrix if
$b_{ij}\leq 0,i\neq j$, and one of the following equivalent conditions holds:

  1. (1)

$B$ is invertible, and $B^{-1}\geq 0$;

  2. (2)

the leading principal minors of the matrix $B$ are positive.

######  Remark 2.2.

Note that for any three column vectors $X\leq Y$ and $Z\geq 0$ we have
$Z^{T}(Y-X)\geq 0$ as a sum of non-negative numbers. Thus, for any matrix
$A\geq 0$ and $X\leq Y$ with entries of arbitrary signs, we have $AX\leq AY$.

For a fixed bounded interval $\Omega=[t_{0},t_{1}]$, let
$\|y\|_{\Omega}=\mbox{ess}\sup_{t\in\Omega}\|y(t)\|$, also for a half-line
$\|f\|_{[t_{0},\infty)}=\mbox{ess}\sup_{t\geq t_{0}}\|f(t)\|$. Denote by
${{\bf L}}_{\infty}^{d}[t_{0},t_{1}]$ the space of all essentially bounded on
$\Omega$ vector functions $y$ with the norm $\|y\|_{\Omega}$, the same for
${{\bf L}}_{\infty}^{d}[t_{0},\infty)$ and for the spaces ${{\bf
L}}_{\infty}^{d\times d}[t_{0},t_{1}]$, ${{\bf L}}_{\infty}^{d\times
d}[t_{0},\infty)$ of matrix functions.

Assume that for a fixed $t_{0}\geq 0$ vector-functions
$\varphi_{i}\colon(-\infty,t_{0}]\rightarrow{\mathbb{R}}^{d}$, $i=1,\dots,n$
are Borel measurable and bounded. Below for every $t_{0}\geq 0$, system (2.1)
will be considered with the initial condition

(2.2) $x_{i}(t)=\varphi_{i}(t),\,\,\,\,t\leq t_{0},~{}~{}i=1,\dots,n$

under the hypotheses:

  * $(i)$

$A_{i},B_{ij}^{k},K_{ij}\colon[0,\infty)\to\mathbb{R}^{d\times d}$,
$k=1,\dots,m_{ij}$, $i,j=1,\dots,n$ belong to ${{\bf L}}_{\infty}^{d\times
d}[0,\infty)$.

  * $(ii)$

$h_{i},h_{ij}^{k},g_{ij}\colon[0,\infty)\to\mathbb{R}$, are Lebesgue
measurable functions, and there are numbers
$\tau_{ij}^{k}>0,\sigma_{ij}>0,\tau_{i}\geq 0$ such that almost everywhere

(2.3) $0\leq t-h_{ij}^{k}(t)\leq\tau_{ij}^{k},~{}0\leq
t-g_{ij}(t)\leq\sigma_{ij},~{}0\leq
t-h_{i}(t)\leq\tau_{i},~{}k=1,\dots,m_{ij},~{}i,j=1,\dots,n,~{}t\geq 0.$

A solution of problem (2.1), (2.2) is understood in the sense of the following
definition.

######  Definition 2.3.

A set of vector-functions
$x_{i}\colon{\mathbb{R}}\rightarrow{\mathbb{R}}^{d},~{}i=1,\dots,n$ is called
a solution of problem (2.1), (2.2) if $x_{i}$ satisfy (2.1) for almost all
$t\in[t_{0},\infty)$ and (2.2) for all $t\in[t_{0},t_{0}-\tilde{\tau}_{i}]$,
where
$\displaystyle\tilde{\tau}_{i}=\max\left\\{\max_{j,k}\tau_{ij}^{k},\max_{j}\sigma_{ij},\tau_{i}\right\\}$,
$i=1,\dots,n$.

Let $x_{i},i=1,\dots,n$ be a solution of problem (2.1), (2.2). Introduce the
matrix function $\Phi(t)=\\{\varphi_{1}(t),\dots,\varphi_{n}(t)\\}$ with
$\varphi_{i}$ as columns.

######  Definition 2.4.

System (2.1) is called uniformly exponentially stable, if there exist positive
constants $H$ and $\nu$, such that any solution $X$ of (2.1), (2.2) satisfies

$\|x_{i}(t)\|\leq
He^{-\nu(t-t_{0})}\sup_{t\in[t_{0},t_{0}-\max_{i}\tilde{\tau}_{i}]}\|\Phi(t)\|,\,\,\,t\geq
t_{0}\geq 0,~{}~{}i=1,\dots,n,$

where $H$ and $\nu$ do not depend on $t_{0}$ and $\Phi$.

Consider a non-homogeneous version of (2.1)

(2.4) $\begin{array}[]{ll}\dot{x_{i}}(t)=&\displaystyle
A_{i}(t)x_{i}(h_{i}(t))+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(t)x_{j}(h_{ij}^{k}(t))\vspace{2mm}\\\
&\displaystyle+\sum_{j=1}^{n}\int\limits_{g_{ij}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds+f_{i}(t),~{}i=1,\dots,n,~{}t\geq
t_{0}\geq 0\end{array}$

and assume that $f_{i}\in{{\bf L}}_{\infty}^{d}[t_{0},\infty)$. The definition
of a solution for (2.4), (2.2) is similar to that of (2.1), (2.2).

A modification of [4, Lemma 2] will further be used.

######  Lemma 2.5.

Let a function $a\colon[t_{0},\infty)\to[0,\infty)$ be Lebesgue measurable and
$\omega\in{{\bf L}}_{\infty}^{1}[t_{0},\infty)$. Then the inequality
$\displaystyle\left|\int_{t_{0}}^{t}e^{-\int_{s}^{t}a(\tau)d\tau}a(s)\omega(s)ds\right|\leq{\mathrm{ess\,sup}}_{t\in[t_{0},t_{1}]}|\omega(t)|,\,\,\,\,t\in[t_{0},t_{1}],$
holds for any $t_{1}>t_{0}$.

Let $x_{i}$ be a solution of non-homogeneous system (2.4), where $f_{i}\in{\bf
L}_{\infty}^{d}[t_{0},\infty)$ satisfy

(2.5) $x_{i}(t)=0,\,\,\,\,t\leq t_{0},~{}~{}i=1,\dots,n.$

The following Bohl-Perron type result is cited from [2, 11].

######  Lemma 2.6.

If, for any $t_{0}\geq 0$ and any $f_{i}\in{{\bf
L}}_{\infty}^{d}[t_{0},\infty)$, $i=1,\dots,n$, all $x_{i}$ in the solution of
problem (2.4), (2.5) belong to ${{\bf L}}_{\infty}^{d}[t_{0},\infty)$,
$i=1,\dots,n$ then system (2.1) is uniformly exponentially stable.

Below, we use the Coppel inequality [7, 21] for a system of ordinary
differential equations

(2.6) $\dot{x}(t)=D(t)x(t),\,\,\,\,\,t\geq 0,$

where the columns of $D$ are in ${{\bf L}}_{\infty}^{d}[0,\infty)$. Let
$\Phi(t,s)$ be the fundamental matrix of (2.6), i.e. a solution of the problem

$\dot{x}(t)=D(t)x(t),\,\,\,\,\,t>s\geq 0,\,\,\,\,\,x(s)=I.$

######  Lemma 2.7 (Coppel inequality).

The fundamental matrix $\Phi(t,s)$ of (2.6) satisfies

(2.7)
$\left\|\Phi(t,s)\right\|\leq\exp\left({\int_{s}^{t}\mu(D(\xi))d\xi}\right),\,\,\,\,\,t>s\geq
0.$

######  Remark 2.8.

By Lemma 2.7, the condition $\mu(D(t))\leq d_{0}<0$, $t\in[t_{0},\infty)$,
$\forall t_{0}\geq 0$ for the matrix measure implies uniform exponential
stability of system (2.6).

The final auxiliary result gives an a priori estimate for the derivative of a
solution. Let $\Omega:=[t_{0},t_{1}]\subset[t_{0},\infty)$, $t_{1}>t_{0}$ and

$C_{ij}(t):=\int_{g_{ij}(t)}^{t}|K_{ij}(t,s)|ds.$

######  Lemma 2.9.

Let $x_{i}=x_{i}(t)$, $i=1,\dots,n$ be a solution of problem (2.4),(2.5),
$\Omega=[t_{0},t_{1}]$. Then the derivatives $\dot{x_{i}}$ satisfy

(2.8)
${\|\dot{x_{i}}\|}_{\Omega}\leq\|A_{i}\|_{[t_{0},\infty)}{\|x_{i}\|}_{\Omega}+\sum_{j=1}^{n}\left(\sum_{k=1}^{m_{ij}}\|B_{ij}^{k}\|_{[t_{0},\infty)}+\|C_{ij}\|_{[t_{0},\infty)}\right){\|x_{j}\|}_{\Omega}+\|f_{i}\|_{[t_{0},\infty)}.$

###### Proof.

Let $t\in\Omega$ and $x=x(t)$ be a solution of (2.4), (2.5). By (2.5),
$x(t)=0$ if $t\leq t_{0}$. Estimating the norm in the left-hand side of (2.4),
we derive

(2.9)
$\|\dot{x_{i}}(t)\|\leq\|A_{i}\|_{[t_{0},\infty)}{\|x_{i}\|}_{\Omega}+\sum_{j=1}^{n}\left(\sum_{k=1}^{m_{ij}}\|B_{ij}^{k}\|_{[t_{0},\infty)}+\|C_{ij}\|_{[t_{0},\infty)}\right){\|x_{j}\|}_{\Omega}+\|f_{i}\|_{\Omega}.$

Since the right-hand side of (2.9) does not depend on $t\in\Omega$, inequality
(2.9) implies (2.8). ∎

##  3\. Main Results

###  3.1. The main theorem

######  Theorem 3.1.

Let there exist $\alpha_{i}$ such that $\mu(A_{i}(t))\leq\alpha_{i}<0$,
$i=1,\dots,n$ and $L=I-D$ be an $M-$matrix, where the entries of
$D=(d_{ij})_{i,j=1}^{n}$ be defined as follows:

$\begin{array}[]{ll}d_{ii}=&\displaystyle\tau_{i}\left\|\frac{A_{i}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\left(\|A_{i}\|_{[t_{0},\infty)}+\sum_{k=1}^{m_{ii}}\|B_{ii}^{k}\|_{[t_{0},\infty)}+\|C_{ii}\|_{[t_{0},\infty)}\right)\displaystyle+\sum_{k=1}^{m_{ii}}\left\|\frac{B_{ii}^{k}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\\!\\!\\!\\!\\!\\!+\left\|\frac{C_{ii}}{\mu(A_{i})}\right\|_{[t_{0},\infty)},\vspace{2mm}\\\
d_{ij}=&\displaystyle\tau_{i}\left\|\frac{A_{i}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\left(\sum_{k=1}^{m_{ij}}\|B_{ij}^{k}\|_{[t_{0},\infty)}+\|C_{ij}\|_{[t_{0},\infty)}\right)+\sum_{k=1}^{m_{ij}}\left\|\frac{B_{ij}^{k}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}+\left\|\frac{C_{ij}}{\mu(A_{i})}\right\|_{[t_{0},\infty)},~{}~{}j\neq
i.\end{array}$

Then system (2.1) is uniformly exponentially stable.

###### Proof.

To apply Lemma 2.6, we explore boundedness of solutions to (2.4),(2.5). First,
transform (2.4) to the form

$\dot{x_{i}}(t)=A_{i}(t)x_{i}(t)-A_{i}(t)\int\limits_{h_{i}(t)}^{t}\dot{x_{i}}(s)ds+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(t)x_{j}(h_{ij}^{k}(t))+\sum_{j=1}^{n}\int\limits_{g_{ij}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds+f_{i}(t).$

Denote by $\Phi_{i}(t,s)$ the fundamental matrix of

(3.1) $\dot{x}(t)=A_{i}(t)x(t).$

Since $\mu(A_{i}(t))\leq\alpha_{i}<0$, (3.1) is uniformly exponentially
stable. Continue transformations of equation (2.4)

$\begin{array}[]{ll}x_{i}(t)=&\displaystyle\int_{t_{0}}^{t}\Phi_{i}(t,s)\left[-A_{i}(s)\int\limits_{h_{i}(s)}^{s}\dot{x_{i}}(\xi)d\xi+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(s)x_{j}(h_{ij}^{k}(s))\right.\\\
&\displaystyle\left.+\sum_{j=1}^{n}\int\limits_{g_{ij}(s)}^{s}K_{ij}(s,\xi)x_{j}(\xi)d\xi\right]ds+\tilde{f_{i}}(t),\end{array}$

where $\tilde{f_{i}}(t)=\int_{t_{0}}^{t}\Phi_{i}(t,s)f_{i}(s)ds\in{{\bf
L}}_{\infty}^{d}[t_{0},\infty)$. Hence by Lemma 2.7,

$\begin{array}[]{lll}\|x_{i}(t)\|&\leq&\displaystyle\int_{t_{0}}^{t}\left.\left.e^{\int_{s}^{t}\mu(A_{i}(\xi))d\xi}\right|\mu(A_{i}(s))\right|\left[\left\|\frac{A_{i}(s)}{\mu(A_{i}(s))}\right\|\int\limits_{h_{i}(s)}^{s}\|\dot{x_{i}}(\xi)\|d\xi\right.\vspace{2mm}\\\
&&+\displaystyle\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}\left\|\frac{B_{ij}^{k}(s)}{\mu(A_{i}(s))}\right\|\|x_{j}(h_{ij}^{k}(s))\|\left.\displaystyle+\sum_{j=1}^{n}\left\|\frac{C_{ij}(s)}{\mu(A_{i}(s))}\right\|\max_{\xi\in[g_{ij}(s),s]}\|x_{j}(\xi)\|\right]ds+\|\tilde{f_{i}}(t)\|\vspace{2mm}\\\
&\leq&\displaystyle\tau_{i}\left\|\frac{A_{i}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\|\dot{x_{i}}\|_{\Omega}+\|\tilde{f_{i}}\|_{[t_{0},\infty)}+\sum_{j=1}^{n}\left(\sum_{k=1}^{m_{ij}}\left\|\frac{B_{ij}^{k}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}+\left\|\frac{C_{ij}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\right)\|x_{j}\|_{\Omega}.\end{array}$

Therefore by Lemma 2.9,

$\begin{array}[]{lll}\|x_{i}\|_{\Omega}&\leq&\displaystyle\tau_{i}\left\|\frac{A_{i}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\left(\|A_{i}\|_{[t_{0},\infty)}{\|x_{i}\|}_{\Omega}+\sum_{j=1}^{n}\left(\sum_{k=1}^{m_{ij}}\|B_{ij}^{k}\|_{[t_{0},\infty)}+\|C_{ij}\|_{[t_{0},\infty)}\right){\|x_{j}\|}_{\Omega}\right)\vspace{2mm}\\\
&&\displaystyle+\sum_{j=1}^{n}\left(\sum_{k=1}^{m_{ij}}\left\|\frac{B_{ij}^{k}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}+\left\|\frac{C_{ij}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\right)\|x_{j}\|_{\Omega}+M_{i}\vspace{2mm}\\\
&=&\displaystyle d_{ii}\|x_{i}\|_{\Omega}+\sum_{j\neq
i}d_{ij}\|x_{j}\|_{\Omega}+M_{i},\end{array}$

where $M_{i}>0$ are constants, not dependent on $t_{1}>t_{0}$.

Introduce the column vectors
$X_{\Omega}=(\|x_{1}\|_{\Omega},\dots,\|x_{n}\|_{\Omega})^{T},M=(\|M_{1}\|,\dots,\|M_{n}\|)^{T}$.
Then $(I-D)X_{\Omega}\leq M$. Since $I-D$ is an $M$-matrix, there exists
$(I-D)^{-1}\geq 0$. Therefore (see Remark 2.2) $X_{\Omega}\leq(I-D)^{-1}M$,
where the right-hand side of the last inequality does not depend on the
interval $\Omega$. Hence the solution of problem (2.4),(2.5) is bounded on
$[t_{0},\infty)$. By Lemma 2.6, system (2.1) is uniformly exponentially
stable. ∎

###  3.2. Corollaries of the main theorem

Consider several particular cases of (2.1).

We start with the system where the first term in each vector equation is non-
delayed

(3.2)
$\dot{x_{i}}(t)=A_{i}(t)x_{i}(t)+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(t)x_{j}(h_{ij}^{k}(t))+\sum_{j=1}^{n}\int\limits_{g_{ij}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds,~{}i=1,\dots,n.$

######  Corollary 3.2.

Assume that $\mu(A_{i}(t))\leq\alpha_{i}<0$, $i=1,\dots,n$ and $L=I-D$ is an
$M-$matrix, where $D=(d_{ij})_{i,j=1}^{n}$ has the entries

$d_{ij}=\sum_{k=1}^{m_{ij}}\left\|\frac{B_{ij}^{k}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}+\left\|\frac{C_{ij}}{\mu(A_{i})}\right\|_{[t_{0},\infty)},~{}i,j=1,\dots,n.$

Then system (3.2) is uniformly exponentially stable.

Next, consider a system without integral terms and with single delays in the
diagonal terms

(3.3) $\dot{x_{i}}(t)=A_{i}(t)x_{i}(h_{i}(t))+\sum_{j\neq
i}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(t)x_{j}(h_{ij}^{k}(t)),~{}i=1,\dots,n.$

######  Corollary 3.3.

Assume that $\mu(A_{i}(t))\leq\alpha_{i}<0,i=1,\dots,n$, and the matrix $I-D$
is an $M-$matrix, where $D=(d_{ij})_{i,j=1}^{n}$ are defined as

$d_{ii}=\displaystyle\tau_{i}\left\|\frac{A_{i}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\|A_{i}\|_{[t_{0},\infty)},~{}d_{ij}=\tau_{i}\left\|\frac{A_{i}}{\mu(A_{i})}\right\|_{[t_{0},\infty)}\sum_{k=1}^{m_{ij}}\|B_{ij}^{k}\|_{[t_{0},\infty)}+\sum_{k=1}^{m_{ij}}\left\|\frac{B_{ij}^{k}}{\mu(A_{i})}\right\|_{[t_{0},\infty)},~{}j\neq
i.$

Then system (3.3) is uniformly exponentially stable.

We proceed to the case of two vector equations:

(3.4)
$\begin{array}[]{l}\displaystyle\dot{x_{1}}(t)=A_{1}(t)x_{1}(h_{1}(t))+\sum_{j=1}^{2}\sum_{k=1}^{m_{1j}}B_{1j}^{k}(t)x_{j}(h_{1j}^{k}(t))+\sum_{j=1}^{2}\int_{g_{1j}(t)}^{t}K_{1j}(t,s)x_{j}(s)ds,\vspace{2mm}\\\
\displaystyle\dot{x_{2}}(t)=A_{2}(t)x_{2}(h_{2}(t))+\sum_{j=1}^{2}\sum_{k=1}^{m_{2j}}B_{2j}^{k}(t)x_{j}(h_{2j}^{k}(t))+\sum_{j=1}^{2}\int_{g_{2j}(t)}^{t}K_{2j}(t,s)x_{j}(s)ds.\end{array}$

######  Corollary 3.4.

Let $\mu(A_{i}(t))\leq\alpha_{i}<0,i=1,2$ and

(3.5) $d_{11}<1,~{}d_{22}<1,~{}d_{12}d_{21}<(1-d_{11})(1-d_{22}),$

where

$\begin{array}[]{ll}d_{11}=&\displaystyle\tau_{1}\left\|\frac{A_{1}}{\mu(A_{1})}\right\|_{[t_{0},\infty)}\left(\|A_{1}\|_{[t_{0},\infty)}+\sum_{k=1}^{m_{11}}\left\|B_{11}^{k}\right\|_{[t_{0},\infty)}\\!\\!\\!\\!\\!\\!+\left\|C_{11}\right\|_{[t_{0},\infty)}\right)+\sum_{k=1}^{m_{11}}\left\|\frac{B_{11}^{k}}{\mu(A_{1})}\right\|_{[t_{0},\infty)}\\!\\!\\!\\!\\!\\!+\left\|\frac{C_{11}}{\mu(A_{1})}\right\|_{[t_{0},\infty)},\\\
d_{22}=&\displaystyle\tau_{2}\left\|\frac{A_{2}}{\mu(A_{2})}\right\|_{[t_{0},\infty)}\left(\|A_{2}\|_{[t_{0},\infty)}+\sum_{k=1}^{m_{22}}\left\|B_{22}^{k}\right\|_{[t_{0},\infty)}\\!\\!\\!\\!\\!\\!+\left\|C_{22}\right\|_{[t_{0},\infty)}\right)+\sum_{k=1}^{m_{22}}\left\|\frac{B_{22}^{k}}{\mu(A_{2})}\right\|_{[t_{0},\infty)}\\!\\!\\!\\!\\!\\!+\left\|\frac{C_{22}}{\mu(A_{2})}\right\|_{[t_{0},\infty)},\end{array}$
$\begin{array}[]{ll}d_{12}=&\displaystyle\tau_{1}\left\|\frac{A_{1}}{\mu(A_{1})}\right\|_{[t_{0},\infty)}\left(\sum_{k=1}^{m_{12}}\left\|B_{12}^{k}\right\|_{[t_{0},\infty)}+\left\|C_{12}\right\|_{[t_{0},\infty)}\right)+\displaystyle\sum_{k=1}^{m_{12}}\left\|\frac{B_{12}^{k}}{\mu(A_{1})}\right\|_{[t_{0},\infty)}+\left\|\frac{C_{12}}{\mu(A_{1})}\right\|_{[t_{0},\infty)},\\\
d_{21}=&\displaystyle\tau_{2}\left\|\frac{A_{2}}{\mu(A_{2})}\right\|_{[t_{0},\infty)}\left(\sum_{k=1}^{m_{21}}\left\|B_{21}^{k}\right\|_{[t_{0},\infty)}+\left\|C_{21}\right\|_{[t_{0},\infty)}\right)+\sum_{k=1}^{m_{21}}\left\|\frac{B_{21}^{k}}{\mu(A_{2})}\right\|_{[t_{0},\infty)}+\left\|\frac{C_{21}}{\mu(A_{1})}\right\|_{[t_{0},\infty)}.\end{array}$

Then system (3.4) is uniformly exponentially stable.

For only two delays in each vector equation in (3.4), we get

(3.6)
$\begin{array}[]{lll}\dot{x_{1}}(t)&=&A_{1}(t)x_{1}(h_{1}(t))+B_{12}(t)x_{2}(h_{12}(t)),\\\
\dot{x_{2}}(t)&=&A_{2}(t)x_{2}(h_{2}(t))+B_{21}(t)x_{1}(h_{21}(t)).\end{array}$

######  Corollary 3.5.

Let $\mu(A_{i}(t))\leq\alpha_{i}<0$, $i=1,2$ and there exist constant matrices
$\tilde{A_{1}},\tilde{A_{2}}$, $\tilde{B}_{12}$, $\tilde{B}_{21}$, such that
$|A_{i}(t)|\leq\tilde{A_{i}}$, $i=1,2$, $|B_{12}(t)|\leq\tilde{B}_{12}$,
$|B_{21}(t)|\leq\tilde{B}_{21}$. If

(3.7)
$\begin{array}[]{l}\displaystyle\tau_{i}\|\tilde{A_{i}}\|^{2}<|\alpha_{i}|,i=1,2,\vspace{2mm}\\\
\displaystyle\|\tilde{B}_{12}\|\|\tilde{B}_{21}\|(1+\tau_{1}\|\tilde{A_{1}}\|)(1+\tau_{2}\|\tilde{A_{2}}\|)<(|\alpha_{1}|-\tau_{1}\|\tilde{A_{1}}\|^{2})(|\alpha_{2}|-\tau_{2}\|\tilde{A_{2}}\|^{2})\end{array}$

then system (3.6) is uniformly exponentially stable.

###### Proof.

Let us check that conditions of Corollary 3.4 hold, where $B_{ii}\equiv
0,C_{ij}\equiv 0$. For brevity, we omit the interval index $[t_{0},\infty)$
for the norms of matrix functions. We have

$d_{ii}=\tau_{i}\left\|\frac{A_{i}}{\mu(A_{i})}\right\|\|A_{i}\|,~{}i=1,2,~{}$
$d_{12}=\tau_{1}\left\|\frac{A_{1}}{\mu(A_{1})}\right\|\|B_{12}\|+\left\|\frac{B_{12}}{\mu(A_{1})}\right\|,~{}d_{21}=\tau_{2}\left\|\frac{A_{2}}{\mu(A_{2})}\right\|\|B_{21}\|+\left\|\frac{B_{21}}{\mu(A_{2})}\right\|.$

The inequality $|A|\leq B$ implies $\||A|\|\leq\|B\|$ for any matrix norm.
Hence

(3.8) $d_{ii}\leq\tau_{i}\frac{\|\tilde{A_{i}}\|^{2}}{|\alpha_{i}|},i=1,2,$
(3.9)
$d_{12}\leq\frac{\|\tilde{B}_{12}\|(1+\tau_{1}\|\tilde{A_{1}}\|)}{|\alpha_{1}|},~{}d_{21}\leq\frac{\|\tilde{B}_{21}\|(1+\tau_{2}\|\tilde{A_{2}}\|)}{|\alpha_{2}|}.$

Inequality (3.8) and the first inequality in (3.7) imply $d_{ii}<1$, $i=1,2$.
By (3.9), inequality $d_{12}d_{21}<(1-d_{11})(1-d_{22})$ holds if

$\frac{\|\tilde{B}_{12}\|\|\tilde{B}_{21}\|(1+\tau_{1}\|\tilde{A_{1}}\|)(1+\tau_{2}\|\tilde{A_{2}}\|)}{|\alpha_{1}||\alpha_{2}|}<\left(1-\frac{\tau_{1}\|\tilde{A_{1}}\|^{2}}{|\alpha_{1}|}\right)\left(1-\frac{\tau_{2}\|\tilde{A_{2}}\|^{2}}{|\alpha_{2}|}\right),$

which is equivalent to the second inequality in (3.7). Since all the
conditions in (3.5) hold, by Corollary 3.4 system (3.6) is uniformly
exponentially stable. ∎

######  Remark 3.6.

Theorem 3.1 and its corollaries also hold for vector differential equations if
the terms $\displaystyle\int\limits_{g_{ij}^{k}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds$
in equation (2.1) are replaced with $\displaystyle
C_{ij}(t)\int\limits_{g_{ij}(t)}^{t}d_{s}R_{ij}(t,s)x_{j}(s)$, where
distributed delays are of a more general type, and
$\displaystyle\int\limits_{g_{ij}(t)}^{t}d_{s}|R_{ij}(t,s)|=I$.

##  4\. Applications to Second Order Vector Delay Differential Equations

###  4.1. Linear equations

In this section we consider a second order vector equation

(4.1) $\ddot{x}(t)+A(t)\dot{x}(t)+B(t)x(h(t))=0,~{}~{}t\geq t_{0}\geq 0,$

where for the matrices $A,B$ and the delay $h$ the same conditions hold as
$(i),(ii)$ for system (2.1). Definitions of solutions of an initial value
problem for vector equation (4.1) and of uniform exponential stability are
also similar to the definitions for (2.1), and we omit them.

######  Theorem 4.1.

Assume that for some $t\geq t_{0}\geq 0$, there exist a constant matrix
$\tilde{A}$ and numbers $\alpha<0$, $\tau\geq 0$ such that $t-h(t)\leq\tau$,
$\mu(-\tilde{A})<0$, $\mu\left(\tilde{A}-2A(t)\right)\leq\alpha<0$ and

(4.2)
$\frac{1}{|\mu(-\tilde{A})|}\left(\left\|\frac{2\tilde{A}A-\tilde{A}^{2}-4B}{\mu(\tilde{A}-2A)}\right\|_{[t_{0},\infty)}\\!\\!\\!\\!\\!\\!+2\tau\left\|\frac{\tilde{A}B}{\mu(\tilde{A}-2A)}\right\|_{[t_{0},\infty)}\right)<1-2\tau\left\|\frac{B}{\mu(\tilde{A}-2A)}\right\|_{[t_{0},\infty)}\\!\\!.$

Then equation (4.1) is uniformly exponentially stable.

###### Proof.

Equation (4.1) can be transformed to the form

(4.3) $\ddot{x}(t)+A(t)\dot{x}(t)+B(t)x(t)-B(t)\int_{h(t)}^{t}\dot{x}(s)ds=0.$

After the substitution

$\dot{x}(t)=-\frac{\tilde{A}}{2}x(t)+y(t),~{}~{}\ddot{x}(t)=\frac{\tilde{A}^{2}}{4}x(t)-\frac{\tilde{A}}{2}y(t)+\dot{y}(t),$

equation (4.3) becomes

$\frac{\tilde{A}^{2}}{4}x(t)-\frac{\tilde{A}}{2}y(t)+\dot{y}(t)+A(t)\left[-\frac{\tilde{A}}{2}x(t)+y(t)\right]+B(t)x(t)-B(t)\int\limits_{h(t)}^{t}\left(-\frac{\tilde{A}}{2}x(s)+y(s)\right)ds=0.$

Hence equation (4.1) is equivalent to the system

(4.4)
$\begin{array}[]{ll}\dot{x}(t)=&\displaystyle-\frac{\tilde{A}}{2}x(t)+y(t)\vspace{2mm}\\\
\dot{y}(t)=&\displaystyle\left(\frac{\tilde{A}}{2}-A(t)\right)y(t)+\left(\frac{\tilde{A}}{2}A(t)-\frac{\tilde{A}^{2}}{4}-B(t)\right)x(t)\vspace{2mm}\\\
&\displaystyle-\frac{\tilde{A}}{2}B(t)\int_{h(t)}^{t}x(s)ds+B(t)\int_{h(t)}^{t}y(s)ds=0.\end{array}$

For system (4.4), let us apply Corollary 3.4. System (4.4) has a form of
(3.4), where

$x_{1}=x,~{}x_{2}=y,~{}~{}A_{1}=-\frac{\tilde{A}}{2},~{}~{}h_{1}(t)=t,~{}\tau_{1}=0,~{}~{}m_{11}=0,~{}m_{12}=1,$
$B_{11}=0,~{}B_{12}=1,~{}~{}h_{12}(t)=t,~{}~{}C_{11}=C_{12}=0,$
$A_{2}=\frac{\tilde{A}}{2}-A(t),~{}h_{2}(t)=t,~{}\tau_{2}=0,~{}~{}m_{21}=1,~{}m_{22}=0,~{}~{}B_{21}=\frac{\tilde{A}}{2}A(t)-\frac{\tilde{A}^{2}}{4}-B(t),$
$C_{21}=\frac{\tilde{A}}{2}|B(t)|(t-h(t))\leq\frac{\tilde{A}}{2}|B(t)|\tau,~{}~{}C_{22}=|B(t)|(t-h(t))\leq|B(t)|\tau.$

We have

$\mu(A_{1})<0,~{}\mu(A_{2})<0,~{}~{}d_{11}=0,~{}d_{22}\leq\tau\left\|\frac{B}{\mu(\frac{\tilde{A}}{2}-A)}\right\|_{[t_{0},\infty)},d_{12}=\frac{2}{|\mu(-\tilde{A})|},$
$d_{21}\leq\left\|\frac{\frac{\tilde{A}}{2}A-\frac{\tilde{A}^{2}}{4}-B}{\mu(\frac{\tilde{A}}{2}-A)}\right\|_{[t_{0},\infty)}+\tau\left\|\frac{\frac{\tilde{A}}{2}B}{\mu(\frac{\tilde{A}}{2}-A)}\right\|_{[t_{0},\infty)}.$

Then $d_{12}d_{21}\leq(1-d_{11})(1-d_{22})$ takes the form

$\frac{2}{|\mu(-\tilde{A})|}\left(\left\|\frac{\frac{\tilde{A}}{2}A-\frac{\tilde{A}^{2}}{4}-B}{\mu(\frac{\tilde{A}}{2}-A)}\right\|_{[t_{0},\infty)}+\tau\left\|\frac{\frac{\tilde{A}}{2}B}{\mu(\frac{\tilde{A}}{2}-A)}\right\|_{[t_{0},\infty)}\right)<1-\tau\left\|\frac{B}{\mu(\frac{\tilde{A}}{2}-A)}\right\|_{[t_{0},\infty)},$

which is equivalent to (4.2). Hence all the conditions of Corollary 3.4 hold
for (4.4). Thus, system (4.4) and therefore equation (4.3) are uniformly
exponentially stable. ∎

######  Remark 4.2.

If $A_{1}\leq A(t)\leq A_{2}$ then anyone of $A_{1},A_{2}$ can be taken as
$\tilde{A}$, see examples below.

######  Corollary 4.3.

Let there exist a constant matrix $\tilde{A}$, $\alpha<0$ and $\tau\geq 0$
such that for $t\geq t_{0}\geq 0$, $t-h(t)\leq\tau$, $\mu(-\tilde{A})<0$,
$\displaystyle\mu\left(\tilde{A}-2A(t)\right)\leq\alpha<0$ and

(4.5)
$\|2\tilde{A}A-\tilde{A}^{2}-4B\|_{[t_{0},\infty)}+2\tau\|\tilde{A}B\|_{[t_{0},\infty)}<|\mu(-\tilde{A})|\left(|\alpha|-2\tau\|B\|_{[t_{0},\infty)}\right).$

Then equation (4.1) is uniformly exponentially stable.

Consider now a second order system with constant coefficients

(4.6) $\ddot{x}(t)+A\dot{x}(t)+Bx(h(t))=0.$

Choosing $\tilde{A}=A$, $\alpha=\mu(-A)$ in Corollary 4.3, we get a result for
(4.6).

######  Corollary 4.4.

Let $t-h(t)\leq\tau<\infty$ for $t\geq t_{0}\geq 0$, $\mu(-A)<0$ and

(4.7) $\|A^{2}-4B\|+2\tau\|AB\|<|\mu(-A)|(|\mu(-A)|-2\tau\|B\|).$

Then equation (4.6) is uniformly exponentially stable.

Next, taking $h(t)=t$ in Theorem 4.1, we get a result for a vector ordinary
differential equation of the second order.

######  Corollary 4.5.

Let there exist a constant matrix $\tilde{A}$, $\alpha<0$ and $\tau\geq 0$
such that for $t\geq t_{0}\geq 0$, $\mu(-\tilde{A})<0$,
$\displaystyle\mu\left(\tilde{A}-2A(t)\right)\leq\alpha<0$ and

(4.8)
$\|2\tilde{A}A-\tilde{A}^{2}-4B\|_{[t_{0},\infty)}<|\alpha||\mu(-\tilde{A})|.$

Then equation (1.1) is uniformly exponentially stable.

######  Corollary 4.6.

Assume that $\mu(-A)<0$ and $\|A^{2}-4B\|<(\mu(-A))^{2}$. Then equation (1.1)
with constant matrices $A$ and $B$ is uniformly exponentially stable.

######  Example 4.7.

Consider equation (4.1), where $t-h(t)\leq\tau$, $\displaystyle
A(t)=\left(\begin{array}[]{cc}4&\sin^{2}t\\\ \cos^{2}t&6\end{array}\right)$,
$\displaystyle B(t)=\left(\begin{array}[]{cc}4&2\sin^{2}t\\\
2\cos^{2}t&8\end{array}\right)$. Choosing
$\displaystyle\tilde{A}:=\left(\begin{array}[]{cc}4&1\\\
1&6\end{array}\right)$, we have $\mu(-\tilde{A})=-3<0$.

$\tilde{A}-2A(t)=\left(\begin{array}[]{cc}-4&1-2\sin^{2}t\\\
1-2\cos^{2}t&-6\end{array}\right),~{}~{}\mu(\tilde{A}-2A(t))\leq\alpha=-3<0$.
Next,

$2\tilde{A}A(t)=\left(\begin{array}[]{cc}32+2\cos^{2}t&12+8\sin^{2}t\\\
8+12\cos^{2}t&72+2\sin^{2}t\end{array}\right),~{}~{}\tilde{A}^{2}=\left(\begin{array}[]{cc}17&10\\\
10&37\end{array}\right),~{}~{}4B(t)=\left(\begin{array}[]{cc}16&8\sin^{2}t\\\
8\cos^{2}t&32\end{array}\right),$
$2\tilde{A}A(t)-\tilde{A}^{2}-4B(t)=\left(\begin{array}[]{ll}-1+2\cos^{2}t&2\\\
-2+4\cos^{2}t&3+2\sin^{2}t\end{array}\right),~{}~{}\tilde{A}B(t)=\left(\begin{array}[]{ll}16+2\cos^{2}t&8+8\sin^{2}t\\\
4+12\cos^{2}t&48+2\sin^{2}t\end{array}\right),$
$\|2\tilde{A}A(t)-\tilde{A}^{2}-4B(t)\|_{[0,\infty)}=7,~{}\|\tilde{A}B\|_{[0,\infty)}=66,\|B\|_{[0,\infty)}=10.$

For $\alpha=-3$, inequality (4.5) is $7+132\tau<3(3-20\tau)$, which holds if
$\tau<\frac{1}{96}$. Thus for $\tau<\frac{1}{96}$ vector equation (4.1) is
uniformly exponentially stable.

Next, consider equation (4.1) with a non-zero right-hand side

(4.9) $\ddot{x}(t)+A(t)\dot{x}(t)+B(t)x(h(t))=L(t),~{}~{}t\geq t_{0}\geq 0.$

The following statement can be justified for a variety of uniformly
exponentially stable differential equations with bounded delays, we prove it
here for second order vector equations.

######  Lemma 4.8.

Let equation (4.1) be uniformly exponentially stable and $L(t)$ have one of
the following properties:

a) $L(t)$ is bounded, $t\in[t_{0},\infty)$;

b) $\int_{t}^{t+1}\|L(s)\|ds$ is bounded, $t\in[t_{0},\infty)$;

c) $\lim_{t\rightarrow\infty}L(t)=0$;

d) $\lim_{t\rightarrow\infty}\int_{t}^{t+1}\|L(s)\|ds=0$.

Then any solution of equation (4.9) possesses the same property.

###### Proof.

According to [5, page 3], [1], there are matrix functions
$X_{1}(t),X_{2}(t),X(t,s)$ such that a solution $x$ of equation (4.9) with the
initial condition $x(t)=\varphi(t)$, $t\leq t_{0}$ satisfies

$x(t)=X_{1}(t)x(t_{0})+X_{2}(t)\dot{x}(t_{0})-\int_{t_{0}}^{t_{0}+\tau}X(t,s)B(s)\varphi(h(s))ds+\int_{t_{0}}^{t}X(t,s)L(s)ds,$

where $\varphi(h(s))=0,~{}h(s)\geq t_{0}$. Exponential stability of (4.1)
implies for some $M_{i}>0$, $\lambda_{i}>0$, $i=1,2$, $M>0$, $\lambda>0$,

$\|X_{i}(t)\|\leq M_{i}e^{-\lambda_{i}t},~{}~{}i=1,2,~{}~{}\|X(t,s)\|\leq
Me^{-\lambda(t-s)}.$

Assume that condition b) holds: $\int_{t}^{t+1}\|L(s)\|ds\leq K$ for $t\geq
t_{0}$. We have

$\begin{array}[]{lll}\|x(t)\|&\leq&\displaystyle
M_{1}e^{-\lambda_{1}t}\|x(t_{0})\|+M_{2}e^{-\lambda_{2}t}\|\dot{x}(t_{0})\|+M\|B\|_{[t_{0},\infty)}\|\varphi\|_{[t_{0}-\tau,t_{0}]}\int\limits_{t_{0}}^{t_{0}+\tau}e^{-\lambda(t-s)}ds\\\
&&\displaystyle+M\int\limits_{t_{0}}^{t}e^{-\lambda(t-s)}\|L(s)\|ds\\\
&\leq&\displaystyle
M_{1}\|x(t_{0})\|+M_{2}\|\dot{x}(t_{0})\|+\frac{M}{\lambda}\|B\|_{[t_{0},\infty)}\|\varphi\|_{[t_{0}-\tau,t_{0}]}e^{\lambda(t_{0}+\tau)}\\\
&&\displaystyle+M\int\limits_{t_{0}}^{t}e^{-\lambda(t-s)}\|L(s)\|ds.\end{array}$

Let $t\in[t_{0}+m,t_{0}+m+1]$ for some integer $m\geq 0$. Then

$\begin{array}[]{ll}\displaystyle\int\limits_{t_{0}}^{t}e^{-\lambda(t-s)}\|L(s)\|ds&\leq\displaystyle\sum_{k=0}^{m}\int\limits_{t_{0}+k}^{t_{0}+k+1}e^{-\lambda(t-s)}\|L(s)\|ds\leq\sum_{k=0}^{m}e^{-\lambda(m-k)}\int\limits_{t_{0}+k}^{t_{0}+k+1}\|L(s)\|ds\\\
&\displaystyle\leq
K\sum_{k=0}^{m}e^{-\lambda(m-k)}<\frac{K}{1-e^{-\lambda}},\end{array}$

where the last expression does not depend on $m$. Hence the solution $x$ of
equation (4.9) is bounded, and $\int_{t}^{t+1}\|x(s)\|\,ds$ is bounded.
Therefore a) and b) yield that solutions of equation (4.9) are bounded.
Similarly, c)-d) imply convergence of all solutions to zero. ∎

######  Theorem 4.9.

Let all the conditions of Theorem 4.1 hold and $L(t)$ have one of the
properties a)-d) of Lemma 4.8. Then any solution of equation (4.9) possesses
the same property.

###  4.2. Nonlinear equations

In this section consider a nonlinear equation

(4.10)
$\ddot{x}(t)+A(t,x(t),\dot{x}(t))\dot{x}(t)+B(t,x(t),\dot{x}(t))x(h(t))=L(t),~{}~{}t\geq
t_{0}\geq 0,$

where the matrices $A(\cdot,u,v)$ and $B(\cdot,u,v)$ belong to ${{\bf
L}}_{\infty}^{d\times d}[0,\infty)$ for any pair $(u,v)$ of locally integrable
on $[0,\infty)$ vector functions, and $L\in{{\bf L}}_{\infty}^{d}[0,\infty)$.
A solution of an initial value problem for (4.10) is defined similarly to that
of (4.1). We assume that a solution exists and is unique.

Consider first equation (4.10) with $L(t)\equiv 0$.

######  Theorem 4.10.

Let for any $t\geq t_{0}\geq 0$, $L(t)\equiv 0$, there exist a constant matrix
$\tilde{A}$ and $\alpha<0$, $\tau\geq 0$ such that $t-h(t)\leq\tau$,
$\mu(-\tilde{A})<0$, and for any pair $(u,v)$ of constant vectors,
$\mu(\tilde{A}-2A(t,u,v))\leq\alpha<0$ and

(4.11)
$\begin{array}[]{ll}&\displaystyle\frac{1}{|\mu(-\tilde{A})|}\left(\left\|\frac{2\tilde{A}A(\cdot,u,v)-\tilde{A}^{2}-4B(\cdot,u,v)}{\mu(\tilde{A}-2A(\cdot,u,v))}\right\|_{[t_{0},\infty)}\right.\vspace{1mm}\\\
+&\left.2\tau\left\|\frac{\tilde{A}B(\cdot,u,v)}{\mu(\tilde{A}-2A(\cdot,u,v))}\right\|_{[t_{0},\infty)}\right)<1-2\tau\left\|\frac{B(\cdot,u,v)}{\mu(\tilde{A}-2A(\cdot,u,v))}\right\|_{[t_{0},\infty)}.\end{array}$

Then any solution $x$ of (4.10) with $L\equiv 0$ tends to zero as
$t\rightarrow\infty$, together with its derivative.

###### Proof.

Let $x$ be a solution of (4.10) with $L\equiv 0$. Consider the linear equation
with variable coefficients for this fixed $x$:

(4.12)
$\ddot{y}(t)+A(t,x(t),\dot{x}(t))\dot{y}(t)+B(t,x(t),\dot{x}(t))y(h(t))=0,~{}~{}t\geq
t_{0}\geq 0.$

All conditions of Theorem 4.1 hold for (4.12). Hence this equation is
uniformly exponentially stable. Then any solution $y$ of this equation tends
to zero, together with its derivative . But the vector-function $x$ is one of
these solutions. Therefore any solution $x$ of (4.10) tends to zero as
$t\rightarrow\infty$, together with its derivative. ∎

Consider now equation (4.10) with a non-zero right-hand side $L(t)$.

######  Corollary 4.11.

Let all the conditions of Theorem 4.10 hold and $L(t)$ possess one of
Properties a)-d) of Lemma 4.8. Then any solution of equation (4.10) has the
same property.

###### Proof.

For any solution $x$ of (4.10), consider the linear equation

(4.13)
$\ddot{y}(t)+A(t,x(t),\dot{x}(t))\dot{y}(t)+B(t,x(t),\dot{x}(t))y(h(t))=L(t),~{}~{}t\geq
t_{0}\geq 0.$

Assume that $L(t)$ satisfies one of a)-d). Since equation (4.13) with $L\equiv
0$ is exponentially stable, by Lemma 4.8 all solutions of (4.13), as well as
the solution $x$ of (4.10), possess this property. ∎

######  Corollary 4.12.

Let there exist a constant matrix $\tilde{A}$, $\alpha<0$ and $\tau\geq 0$
such that $\mu(-\tilde{A})<0$, $t-h(t)\leq\tau$ for $t\geq t_{0}\geq 0$, and
for any pair $(u,v)$ of constant vectors,
$\mu\left(\tilde{A}-2A(t,u,v)\right)\leq\alpha<0$ and

(4.14)
$\begin{array}[]{ll}&\displaystyle\|2\tilde{A}A(\cdot,u,v)-\tilde{A}^{2}-4B(\cdot,u,v)\|_{[t_{0},\infty)}+2\tau\|\tilde{A}B(\cdot,u,v)\|_{[t_{0},\infty)}\\\
<&\displaystyle|\mu(-\tilde{A})|(|\alpha|-2\tau\|B(\cdot,u,v)\|_{[t_{0},\infty)}).\end{array}$

Then any solution $x$ of (4.10) with $L\equiv 0$ tend to zero as
$t\rightarrow\infty$, together with its derivative.

######  Example 4.13.

Consider equation (4.10), where $L(t)\equiv 0$, $t-h(t)\leq\tau$,
$x=(x_{1},x_{2})$,

$A(t)=\left(\begin{array}[]{cc}4&\sin^{2}(tx_{1}(t))\\\
\cos^{2}(tx_{2}(t))&6\end{array}\right),~{}~{}B(t)=\left(\begin{array}[]{cc}4&2\sin^{2}(tx_{1}(t))\\\
2\cos^{2}(tx_{2}(t))&8\end{array}\right).$

Exactly the same calculations as in Example 4.7 and condition (4.14) imply
that for $\tau<\frac{1}{96}$, any solution of (4.10) tends to zero as
$t\rightarrow\infty$, together with its derivative.

##  5\. Discussion and Topics for Further Research

We investigated uniform exponential stability for systems of linear vector
functional differential equations, as well as asymptotic stability of linear
and certain nonlinear vector equations of the second order. System (2.1) has
many applications to real-word models but, to the best of our knowledge, very
little is known about its exponential stability. One of most common
applications is to second order vector delay differential equations, which is
considered in the paper for equations with one delay. However, since the
general theory is developed for a more general model, the results are easily
extended to equations with several bounded delays

(5.1)
$\ddot{x}(t)+A(t)\dot{x}(t)+\sum_{k=1}^{m}B_{k}(t)x(h_{k}(t))=0,~{}~{}t\geq
t_{0}\geq 0,$

and with a bounded distributed delay

(5.2) $\ddot{x}(t)+A(t)\dot{x}(t)+\int_{h(t)}^{t}B(t,s)x(s)ds=0,~{}~{}t\geq
t_{0}\geq 0.$

Corollary 4.5 gives a new simple and general stability test for linear
ordinary differential vector equations of the second order, compared to the
papers [9, 10, 17, 19, 22] where only matrices of a specific form were
considered. For example, in the recent paper [10] the following stability test
was obtained. Denote $\displaystyle
A_{R}=\frac{A+A^{*}}{2},~{}A_{I}=\frac{A-A^{*}}{2i}$ .

######  Proposition 5.1.

Let there exist $m_{A}>0$ such that

$A_{R}(t)\geq 2m_{A}I,~{}T_{R}(t)\geq
m_{A}I,~{}\|T_{R}\|_{[0,\infty)}+\|T_{I}\|_{[0,\infty)}\leq 2m_{A}^{2},$

where $T(t)=m_{A}A(t)-B(t)$ and inequality $A\geq B$ means that $A-B$ is a
positive definite matrix (all eigenvalues are non-negative). Then equation
(1.1) is exponentially stable.

Actually it is not easy to check that all eigenvalues of time-variable matrix
$T(t)$ are non-negative.

For a nonlinear delay vector equation of the second order

(5.3) $\ddot{X}(t)+F(X(t),\dot{X}(t))\dot{X}(t)+H(X(t-\tau))=0,$

where the matrix-function $F(u,v)$ and vector-function $H(u)$ are continuous,
an interesting explicit stability test was obtained in [20].

######  Proposition 5.2.

[20] Let there exist positive numbers $a_{0},a_{1},a_{2}$ such that
1) the matrix $F(u,v)$ is symmetric, and the eigenvalues of this matrix
$\lambda_{i}(F(u,v))\geq a_{1}$ for all pairs
$(u,v)\in{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}$ .
2) $H(0)=0$, $H(X)\neq 0$, $X\neq 0$, the Jacobian $J_{H}(X)$ is symmetric and
$a_{2}\leq\lambda_{i}(J_{H}(X))\leq a_{0}$.

If $\tau<\frac{a_{1}}{a_{0}\sqrt{d}}$ then the trivial solution of (5.3) is
asymptotically stable.

Compare now Proposition 5.2 of [20] and the test in Theorem 4.10.

The stability result of Proposition 5.2 applies only to autonomous equations,
with continuous functions $G,H$ and constant delay $\tau$. Our result works
for non-autonomous equations, with measurable coefficients and a variable
delay. We also apply different stability methods: in [20] the method of
Krasovskii-Lyapunov functionals was applied, in this paper the stability
results are based on Bohl-Perron theorem.

The linear equation

(5.4) $\ddot{X}(t)+A\dot{X}(t)+Bx(t-\tau)=0,$

where $A,B$ are constant matrices, is a partial case of equation (5.3). Let us
compare for (5.4) the results of Corollary 4.4 and Proposition 5.2. For the
scalar case ($d=1$) Proposition 5.2 leads to $B>0$, $A>0$, $B\tau<A$, while
Corollary 4.4 reduces to $B>0,A>0,|A^{2}-4B|<A^{2}-4AB\tau$. In the scalar
case (and presumably also for small $d$), Proposition 5.2 is better than
Corollary 4.4. But for large $d$, Corollary 4.4 could be better than
Proposition 5.2. Assume, for example, that in (5.4) $A=aI,B=bI,a>0,b>0$. Then
the condition of Proposition 5.2 becomes $b\tau<\frac{a}{\sqrt{d}}$, where
$\tau\rightarrow 0$ for fixed $a,b$ and $d\rightarrow\infty$. Proposition 5.2
gives the same condition as for the scalar case: $|a^{2}-4b|<a^{2}-4ab\tau$
and does not depend on $d$. Consider for comparison $A=aI$,
$B=\frac{1}{4}a^{2}I$, where $a>0$. Then Proposition 5.2 implies asymptotic
stability for $\tau<\frac{4}{a\sqrt{d}}$ compared to $\tau<\frac{1}{a}$ in
Corollary 4.4. Thus Corollary 4.4 gives a sharper result for $d>16$.

Together with stability results obtained here, the aim of the paper was also
to attract attention to the new class of systems of functional differential
equations.

It is not yet clear how to obtain exponential stability tests for equations
with unbounded delays, for example, pantograph-type with $h(t)=\lambda t$, and
for equations with delays in the derivative terms such as

(5.5) $\ddot{x}(t)+A(t)\dot{x}(g(t))+B(t)x(h(t))=0.$

Further, let us discuss other possible problems for future research.

An interesting extension of (2.1) is a nonlinear system

$\dot{x_{i}}(t)=A_{i}(t)x_{i}(h_{i}(t))+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}H_{ij}(t,x_{j}(h_{ij}^{k}(t)))+\sum_{j=1}^{n}F_{ij}\left(t,\int\limits_{g_{ij}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds\right),~{}i=1,\dots,n.$

Results for this system can be applied to a nonlinear vector delay equation of
the second order

$\ddot{x}(t)+F(t,\dot{x}(t))+H(t,x(h(t))=0.$

Also, there are no results on asymptotic behavior of solutions to neutral
systems of vector differential equations

$\dot{x_{i}}(t)-Q_{i}(t)\dot{x}_{i}(g_{i}(t))=A_{i}(t)x_{i}(h_{i}(t))+\sum_{j=1}^{n}\sum_{k=1}^{m_{ij}}B_{ij}^{k}(t)x_{j}(h_{ij}^{k}(t))+\sum_{j=1}^{n}\int\limits_{g_{ij}(t)}^{t}K_{ij}(t,s)x_{j}(s)ds,$

$i=1,\dots,n$ and neutral vector equations of the second order

$\ddot{x}(t)+A\dot{x}(g(t))+Bx(h(t))=C\ddot{x}(p(t)).$

It would be interesting to study stability for systems of second order
equations even in the scalar case, such as

$\ddot{x}(t)+A_{1}\dot{x}(g_{1}(t))+B_{1}y(h_{1}(t))=0,$
$\ddot{y}(t)+A_{2}\dot{y}(g_{2}(t))+B_{2}x(h_{2}(t))=0$

and

$\ddot{x}(t)+A_{1}\dot{x}(g_{1}(t))+B_{1}x(h_{1}(t))=C_{1}y(p_{1}(t)),$
$\ddot{y}(t)+A_{2}\dot{y}(g_{2}(t))+B_{2}y(h_{2}(t))=C_{2}x(p_{2}(t)).$

## References

  * [1] N. V. Azbelev, V. Maksimov, L. Rakhmatullina, Introduction to the Theory of Linear Functional-Differential Equations, Advanced Series in Mathematical Science and Engineering,3, World Federation Publishers Company, Atlanta, GA, 1995.
  * [2] N. V. Azbelev, P. M. Simonov, Stability of Differential Equations with Aftereffect, Stability and Control: Theory, Methods and Applications 20, Taylor $\&$ Francis, London, 2003\.
  * [3] L. Berezansky, E. Braverman, On stability of linear neutral differential equations with variable delays. _Czechoslovak Math. J._ 69(144) (2019), 863–891.
  * [4] L. Berezansky, J. Diblík, Z. Svoboda, Z. Šmarda, Simple uniform exponential stability conditions for a system of linear delay differential equations, _Appl. Math. Comput._ 250 (2015), 605–614.
  * [5] L. Berezansky, A. Domoshnitsky, R. Koplatadze, Oscillation, Nonoscillation, Stability and Asymptotic Properties for Second and Higher Order Functional Differential Equations,CRC Press, 2020.
  * [6] A. Berman, R. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Computer Science and Applied Mathematics, Academic Press, New York-London, 1979.
  * [7] C. A. Desoer, M. Vidyasagar, Feedback Systems: Input-Output Properties, Electrical Science, Academic Press, New York, 1975.
  * [8] J. Hale, Theory of Functional Differential Equations. Second edition. Applied Mathematical Sciences, Vol. 3. Springer-Verlag, New York-Heidelberg, 1977.
  * [9] C. J. Harris, J. F. Miles, Stability of Linear Systems: Some Aspects of Kinematic Similarity. Mathematics in Science and Engineering, 153. Academic Press, London-New-York, 1980.
  * [10] M. I. Gil, Stability of linear systems governed by second order vector differential equations, Internat. J. Control 78 (2005), 534–536.
  * [11] M. Gil’, Stability of Vector Differential Delay Equations, Frontiers in Mathematics, Birkhäuser $\&$ Springer Basel AG, Basel, 2013.
  * [12] K. Gopalsamy, Stability and Oscillation in Delay Differential Equations of Population Dynamics, Kluwer Academic Publishers, Dordrecht, Boston, London, 1992.
  * [13] V. B. Kolmanovskii, V. R. Nosov, Stability of functional-differential equations, Mathematics in Science and Engineering, 180. Academic Press, Inc. 1986.
  * [14] V. Kolmanovskii, A. Myshkis, Introduction to the Theory and Applications of Functional-Differential Equations. Mathematics and its Applications, 463. Kluwer Academic Publishers, Dordrecht, 1999.
  * [15] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Mathematics in Science and Engineering 191, Academic Press, Inc., Boston, MA, 1993.
  * [16] P. H. A. Ngoc, Novel criteria for exponential stability of functional differential equations, Proc. Amer. Math. Soc. 141 (2013), no. 9, 3083–3091.
  * [17] W. J. Rugh, Linear System Theory, Prentice Hall Information and System Sciences Series. Prentice Hall, Inc., Englewood Cliffs, NJ, 1993.
  * [18] G. Söderlind, The logarithmic norm. History and modern theory, _BIT_ 46 (2006), no. 3, 631-–652.
  * [19] J. Sun, Q.-G. Wang, Q.-C. Zhong, A less conservative stability test for second-order linear time-varying vector differential equations, Internat. J. Control 80 (2007), 523–526.
  * [20] C. Tunç, Stability to vector Liénard equation with constant deviating argument, Nonlinear Dynam. 73 (2013), 1245–-1251.
  * [21] Z. Zahreddine, Matrix measure and application to stability of matrices and interval dynamical systems, _Int. J. Math. Math. Sci._ 2003:2 (2003), 75–85.
  * [22] A. Zevin, M. Pinsky, Exponential stability and solution bounds for systems with bounded nonlinearities. New directions on nonlinear control, IEEE Trans. Automat. Control 48 (2003), 1799–1804.