# Kinetic instability in inductively oscillatory plasma equilibrium

F. Cruz<EMAIL_ADDRESS>T. Grismayer  L. O. Silva
GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico,
Universidade de Lisboa, 1049-001 Lisboa, Portugal

###### Abstract

A uniform in space, oscillatory in time plasma equilibrium sustained by a
time-dependent current density is analytically and numerically studied
resorting to particle-in-cell simulations. The dispersion relation is derived
from the Vlasov equation for oscillating equilibrium distribution functions,
and used to demonstrate that the plasma has an infinite number of unstable
kinetic modes. This instability represents a new kinetic mechanism for the
decay of the initial mode of infinite wavelength (or equivalently null
wavenumber), for which no classical wave breaking or Landau damping exists.
The relativistic generalization of the instability is discussed. In this
regime, the growth rate of the fastest growing unstable modes scales with
$\gamma_{T}^{-1/2}$, where $\gamma_{T}$ is the largest Lorentz factor of the
plasma distribution. This result hints that this instability is not as
severely suppressed for large Lorentz factor flows as purely streaming
instabilities. The relevance of this instability in inductive electric field
oscillations driven in pulsar magnetospheres is discussed.

††preprint: APS/123-QED

Plasma equilibria are usually found as solutions to a combination of kinetic
or fluid equations and Maxwell’s equations in the stationary limit.
Remarkably, time-dependent equilibrium conditions can also be determined for
certain systems consisting of a plasma and one or more waves. Such systems are
often unstable to parametric instabilities [1], determined by the properties
of the waves supporting the equilibrium. Parametric mode excitation can be
generally understood as the excitation of two or more plasma waves from a pump
wave of finite amplitude $E_{0}$, frequency $\omega_{0}$ and wavenumber
$k_{0}$. The pump can be purely electromagnetic, e.g. a laser [2, 3, 4, 5, 6,
7, 8, 9], of Alfvènic nature [10, 11, 12, 13], or even electrostatic [14, 15].
The scope of application of parametric decay is vast, as unstable parametric
modes have been explored in the context of laser-plasma interactions [7, 8],
inertial [2, 3] and magnetic [9] confinement fusion, laser beam amplification
schemes [16] as well as of space and astrophysical plasmas [17, 18, 10, 19].

In this work, we address the stability of a uniform in space, oscillatory in
time plasma equilibrium. In this equilibrium, $(\omega_{0},k_{0}=0)$ electric
field oscillations, are inductively supported by repeated reversals of the
plasma current. This configuration is similar to the oscillating two-stream
instability 111We refer to the oscillating two-stream instability as the
instability that develops when two uniform plasma populations (namely
electrons and/or positrons) counter-stream with time-dependent average
velocities. This is a configuration similar to that in Ref. [1], but in which
the pump and all other coupled modes have a similar frequency. This should not
be confused with other misleading definitions of the oscillating two-stream
instability, discussed in Ref. [15]., that has been historically addressed as
a parametric instability, but with a time-dependent relative drift velocity
[21]. This is a regime of interest in pulsar magnetospheres, where inductive
oscillations are excited [22, 23] following electron-positron pair cascades in
strong fields [24, 25]. In this regime, a fundamental analytical description
of the instability is difficult to obtain because different $\omega_{0}$
harmonics are coupled. Here, we present a theoretical analysis of this
instability, and show that it acts as a fundamental plasma process for the
transfer of energy from the inductive pump wave to smaller and smaller plasma
kinetic scales.

This Letter is organized as follows: first, we describe the plasma equilibrium
supporting these oscillations. We then present the dispersion relation of
electrostatic plasma waves developed in this equilibrium for an initial
waterbag distribution function, and analyse it both theoretically and
numerically resorting to particle-in-cell (PIC) simulations. Finally, the
conclusions of this work are outlined, and their relevance in pulsar
magnetospheres is discussed.

We consider a uniform unmagnetized pair plasma 222The inductive waves and the
instability presented in this work also develop for a single oscillating
electron population in an ion background. We choose to present the derivation
for pair plasmas due to their applicability to pulsar magnetospheres. in the
presence of a uniform electric field
$\mathbf{E}=E\mathbf{\hat{x}}=E_{0}\mathbf{\hat{x}}$ which we assume, without
loss of generality, to be positive. In the presence of this field, electrons
and positrons are accelerated in opposite directions, driving a current
$\mathbf{j}=j\mathbf{\hat{x}}$ that inductively reduces $E$. The plasma
current is maximum when $E$ vanishes, thus reversing the electric field.
Electrons and positrons are decelerated and $j$ decreases in magnitude until
it is reversed. The inverse process occurs and the system re-establishes the
initial conditions. Since all dynamics is one-dimensional in space, we
hereafter restrict our analysis to the $\mathbf{\hat{x}}$ components of
fields, currents and particle trajectories. To determine the evolution of the
$E$, we first take the time derivative of Ampère’s law,

$\frac{\partial^{2}E}{\partial t^{2}}=-4\pi e\frac{\partial j}{\partial
t}=-8\pi e\int\mathrm{d}p\ v\frac{\partial f_{0}^{+}}{\partial t}\ ,$ (1)

where we have assumed that the plasma current is driven by counter-propagating
positrons and electrons with uniform density $n_{0}$ and average velocity
$\pm\langle v_{+}\rangle$ respectively, i.e. $j=2en_{0}\langle v_{+}\rangle$
($e$ is the elementary charge). We have also used the definition of average
velocity $\langle v_{+}\rangle=\int\mathrm{d}p\ vf_{0}^{+}/\int\mathrm{d}p\
f_{0}^{+}$, where $f_{0}^{+}=f_{0}^{+}(p,t)$ is the positron momentum
distribution function, normalized as $\int\mathrm{d}p\ f_{0}^{+}(p,t)=n_{0}$
(the same applies for the electron distribution function $f_{0}^{-}$). From
the Vlasov equation describing this equilibrium, we can write $\partial
f_{0}^{+}/\partial t=-eE\partial f_{0}^{+}/\partial p$, and perform the
integral in Eq. (1) by parts to obtain

$\frac{\partial^{2}E}{\partial t^{2}}=-\frac{8\pi
e^{2}}{m_{e}}E\int\mathrm{d}p\frac{f_{0}^{+}(t,p)}{\gamma^{3}}\equiv-\bigg{\langle}\frac{8\pi
e^{2}n_{0}}{m_{e}\gamma^{3}}\bigg{\rangle}E\ ,$ (2)

where we have used the relationship between momentum and velocity, $p=\gamma
m_{e}v$, where $m_{e}$ is the electron mass and
$\gamma=1/\sqrt{1-v^{2}/c^{2}}$. The natural oscillation frequency of the
electric field is then $\omega_{0}=\sqrt{\langle 8\pi
e^{2}n_{0}/m_{e}\gamma^{3}\rangle}$. In the non relativistic limit,
$\omega_{0}$ is a constant in time and $E$ is purely harmonic. If the electric
field amplitude is large enough to accelerate particles to relativistic
velocities, the oscillation may be, in general, more complex, with the field
changing purely linearly in time between crests and troughs (i.e. with a
triangular shape) [22]. For simplicity, we present here an analytical
description of these waves in non relativistic regime, and then discuss the
generalization to the relativistic regime. In the non relativistic limit, the
momentum conservation equation of positrons/electrons in the equilibrium
defined in Eq. (2) can be integrated, and their unperturbed orbits are

$v_{\pm}(t)=v_{0}\pm\delta v\cos(\omega_{0}t+\theta)\equiv v_{0}\pm\delta
v\cos\phi\ ,$ (3)

where $\delta v=eE_{0}/m_{e}\omega_{0}$ and $\theta$ is a phase factor such
that $v_{\pm}=v_{0}$ at a reference time $t=t_{0}$. We may look for unstable
plasma modes $\mathbf{k}=k\mathbf{e}_{x}$, with $k\neq 0$, that can develop
and grow exponentially, by integrating the linearized Vlasov equation along
the unperturbed orbits in Eq. (3). This is an approach similar to the
derivation of Bernstein waves [27], extensively documented in the literature
[28, 29]. In this Letter, we discuss the final dispersion relation, and
present all the details of the derivation in the Supplemental Material.

We consider a plasma where the equilibrium distribution function of both
electrons and positrons is a waterbag, $f_{0}^{\pm}(v,t)=n_{0}/\Delta
v(H(v+v_{T}\mp\delta v\cos\phi))-H(v-v_{T}\mp\delta v\cos\phi))$, where
$\Delta v=2v_{T}$ and $H(x)$ is the Heaviside function. The dispersion
relation is

$1-\sum_{n=-\infty}^{+\infty}J_{n}^{2}\left(\frac{k\delta
v}{\omega_{0}}\right)\frac{\omega_{0}^{2}}{(\omega-n\omega_{0})^{2}-k^{2}v_{T}^{2}}=0\
,$ (4)

where $J_{n}(x)$ are Bessel functions of the first kind and order $n$. The
dispersion relation in Eq. (4) readily indicates that an infinite number of
branches $\omega(k)$ exists. These branches correspond to regions in the
$(\omega,k)$ space where each term (in $n$) on the right-hand side of Eq. (4)
dominates the series.

Given the complex form of Eq. (4), general analytical calculations of the
unstable modes and their growth rates are difficult to obtain. Solving Eq. (4)
with each series term individually yields purely real branches,
$\omega_{n\pm}=n\omega_{0}\pm\sqrt{\omega_{0}^{2}J_{n}^{2}+k^{2}v_{T}^{2}}$,
where $J_{n}\equiv J_{n}(k\delta v/\omega_{0})$. The branches with $|n|\leq
2$, relevant for $kv_{T}/\omega_{0}\lesssim 1$ at frequencies near small
multiples of $\omega_{0}$, are plotted in Fig. 1(a) as a function of $k$. An
infinite number of crossings between branches exists, with those between
consecutive branches $(n,n+1)$ occurring at lower wavenumbers. Symmetric
branches $(n,-n)$ cross roughly at $k\simeq n\omega_{0}/v_{T}$ and
$\omega_{r}\equiv\Re(\omega)=0$ (as expected from symmetry). Writing
$\omega=\omega_{r}+i\Gamma$ and assuming that $\Gamma\ll\omega_{0}$, we can
show that the terms of the series in Eq. (4) decrease with $n$, and thus we
can keep only the small $n$ terms of the series to solve the dispersion
relation. Here, we present two analytical solutions of Eq. (4) yielding
unstable modes, corresponding to the interaction between branches i) $n=0,\pm
1$ at small $k$, and ii) $n=\pm 1$ at $kv_{T}/\omega_{0}\simeq 1$. The latter
mode only exists for finite $v_{T}$, whereas the former exists even when
$v_{T}=0$. For this reason, we hereafter refer to these modes as thermal and
fluid, respectively.

To determine the properties of the thermal mode, we look for solutions to Eq.
(4) with $\omega_{r}=0$ with only terms $\pm n$, and then take the particular
case $n=1$. Taking advantage of the Bessel functions’ symmetry
$J_{-n}=(-1)^{n}J_{n}$, we obtain

$\omega^{2}_{n}=n^{2}\omega_{0}^{2}+F(k)\pm\sqrt{\omega_{0}^{4}J_{n}^{4}+4n^{2}\omega_{0}^{2}F(k)}\
,$ (5)

where $F(k)=\omega_{0}^{2}J_{n}^{2}+k^{2}v_{T}^{2}$. The solutions in Eq. (5)
are unstable ($\omega_{n}^{2}<0$) if
$2\omega_{0}^{2}J_{n}^{2}>n^{2}\omega_{0}^{2}-k^{2}v_{T}^{2}$, which is
satisfied for wavenumbers $kv_{T}/\omega_{0}\in(\sqrt{n^{2}-J_{n}^{2}},n)$.
For $n=1$ in particular, the symmetric branches cross at
$kv_{T}/\omega_{0}=1$, and
$\Gamma_{1,\textrm{max}}\equiv\max(\Im(\omega_{1}))\simeq\omega_{0}J_{1}^{2}/2$.
If $k\delta v/\omega_{0}>1$, the large argument asymptotic expansion of Bessel
functions applies, and we find $\Gamma_{1,\textrm{max}}/\omega_{0}\simeq
m_{e}\omega_{0}v_{T}/\pi eE_{0}\simeq v_{T}/\pi\delta v$. In Fig. 1(b), we
plot the imaginary component of $\omega_{n}$ as a function of $k$ for $n=1,2$,
showing that unstable modes indeed occur close to crossings between the
corresponding branches, and that their growth rate decays with $n$.

The fluid mode couples branches $n=0,\pm 1$ at small $k$. For $v_{T}=0$ and
$k\delta v/\omega_{0}\ll 1$, the dispersion relation reduces to

$\frac{(\Omega-1)^{3}}{3\Omega-1}=K\ ,$ (6)

where $\Omega=(\omega/\omega_{0})^{2}$ and $K=(1/2)(k\delta
v/\omega_{0})^{2}$. Eq. (6) has unstable solutions with
$\omega_{r}/\omega_{0}\simeq\pm 1\mp(1/4)(k\delta v/\omega_{0})^{2/3}$ and
growth rate $\Gamma/\omega_{0}\simeq(\sqrt{3}/4)(k\delta v/\omega_{0})^{2/3}$.

Figure 1:  Theoretical prediction of real and imaginary components of wave
frequencies excited with $\delta v/c=0.14$ and $\Delta v/c=0.1$: (a) shows the
purely real solutions obtained from individual branches of the full dispersion
relation in Eq. (4). Lines labeled with $n^{\pm}$ correspond to solutions
where only term $n$ was kept in the series and where
$\partial\omega_{r}/\partial k\gtrless 0$ for $k\gg 1$, respectively; (b)
illustrates the growth rate of unstable modes obtained by combining the
symmetric branches of Eq. (4).

To confirm our theoretical findings, we have performed a set of 1D PIC
simulations with OSIRIS [30, 31] considering a uniform pair plasma of density
$n_{0}$ and a waterbag distribution function in momentum with $\Delta
v/c=0.05-0.3$. The plasma is subject to an initial electric field
$E_{0}/(m_{e}c\omega_{p}/e)=0.2$, where $\omega_{p}^{2}=4\pi
e^{2}n_{0}/m_{e}=\omega_{0}^{2}/2$. The simulation domain has a length
$L/(c/\omega_{0})\simeq 70$, and is discretized in $N=5000$ cells, with $500$
particles/cell/species. The simulation time step is $\Delta
t\omega_{p}=0.005$. When the simulations start, the plasma undergoes the
oscillations described by Eq. (2). Unstable modes grow on top of the $k=0$
oscillations. These oscillations remain initially stable, but are then damped
as energy is transferred to unstable modes and into particle kinetic energy.

Fig. 2 shows the electric field profile and the electron phase space at (a)
the beginning of the simulation, (b) during the linear stage of the
instability and (c) at a time when the instability has saturated for a
simulation with $\Delta v/c=0.1$. In Fig. 2(b), we observe that small
perturbations start growing on top of the oscillating electric field,
modifying the initial waterbag velocity distribution. The initial
electromagnetic energy density ${E_{0}}^{2}/8\pi$ is converted into unstable
modes until they saturate. The initial particle distribution is then strongly
distorted, extending well beyond the initial thermal spread $v_{T}/c=0.05$ for
long times, see Fig. 2(c).

Figure 2:  Temporal evolution of electron phase space (in color) and electric
field (in black lines) for a simulation with
$E_{0}/(m_{e}c\omega_{0}/e)=\delta v/c=0.14$ and $\Delta v/c=0.1$: (a) initial
plasma configuration, (b) linear stage of the instability, showing electric
field perturbations grown on top of its uniform oscillatory component, (c)
final plasma state, after the instability has saturated.

The time evolution of the electric field Fourier spectrum is presented in Fig.
3(a), showing that well-defined unstable modes grow exponentially during the
linear stage of the instability. We observe multiple unstable thermal modes:
the mode with lowest $k$ is the first to grow (region R2), followed by higher
$k$ modes (regions R3-4). We conjecture that these unstable modes are coupled,
as discussed in the Supplemental Material. The overall growth rate and total
energy stored in unstable modes is dominated by the lowest $k$ thermal mode,
as shown in Fig. 3(b), illustrating the electric field energy stored in
regions R1-4 of the $k$ space identified in Fig. 3(a). It is also possible to
observe in Fig. 3(b) that the energy in R1, corresponding to $k\sim 0$,
decreases with time at the expense of the growth of all the other modes. Both
the bandwidth and maximum growth rate of the most unstable modes are in good
agreement with the solutions of Eq. (5). We observe fluid modes grow at
$k\delta v/\omega_{0}\ll 1$ but only weakly at early times, saturating at
levels that do not play any dynamic role on the evolution of the system. When
the unstable thermal modes reach a finite amplitude, particle acceleration
occurs and causes a saturation of the instability. This is followed by
electrostatic turbulence, which causes a strong distortion of the distribution
(see Fig. 2(c)). The numerical dispersion relation of the plasma in this
simulation is presented in Fig. 3(c), which was obtained by Fourier
transforming the data in Fig. 3(a) in time during the linear stage of the
instability. The numerical dispersion relation shows that the inductive mode
is present at $\omega=\omega_{0}$ for low $k$ and that, in general, unstable
modes have a real frequency component multiple of $\omega_{0}$.

Figure 3:  Fourier analysis of unstable modes in the simulation illustrated in
Fig. 2: (a) and (b) show respectively the time evolution of the electric field
Fourier transform and of the energy in $k$ bands corresponding to the main
oscillatory mode (R1) and to the three unstable modes with lowest $k$ (R2-4);
(c) illustrates the numerical plasma dispersion relation, obtained by taking
the Fourier transform of (a) in time during the linear stage of the
instability.

For all simulations performed with $\delta v/v_{T}\gtrsim 1$, the growth rate
of the fastest growing modes is $\Gamma_{\textrm{max}}/\omega_{0}\sim 0.1$ and
decreases with increasing $\delta v/v_{T}$. For $\delta v/v_{T}\lesssim 1$, we
have verified that the growth rate decreases with increasing $v_{T}$. This is
also verified for pair plasmas with Maxwellian distributions with thermal
velocities $v_{\textrm{th}}/c=0.05-0.2$. A detailed discussion of the scaling
of the growth rate with $\delta v/v_{T}$ is presented in the Supplemental
Material. We have found that the fastest growing modes in all simulations with
$\delta v/v_{T}\gtrsim 1$ have $kv_{T}/\omega_{0}\simeq 0.5$, slightly lower
than that predicted from linear theory, $k/(\omega_{0}/v_{T})\simeq 0.6-1$. We
attribute this difference to i) an average performed in the derivation to
obtain a dispersion relation independent of time (see Supplemental Material)
and to ii) the weakly nonlinear regime in which the thermal modes develop.
Regarding ii), we note that fluid modes develop from early times, with which
the distribution may interact and broaden slightly, modifying the properties
of the thermal mode. For $\delta v/v_{T}<1$, simulations show that the fastest
growing modes have $kv_{T}/\omega_{0}\sim\delta v/v_{T}$, a result that the
dispersion relation presented in this work fails to explain. Surprisingly,
however, the growth rate of these modes is well described by linear theory.

We now discuss the generalization to relativistic conditions, i.e.
$E_{0}/(m_{e}c\omega_{0}/e)\gg 1$. In general, under these conditions, the
unperturbed orbits in velocity space become nearly perfect square waves,
oscillating between $\pm c$, i.e. the trajectories $x(t)$ are triangular
waves. It is well known that these can be described as a series of sinusoidal
functions that is well approximated by the lowest order term. The same
approximation can be applied to the electric field profile in regimes where
$\delta p/p_{T}\gg 1$, where $\delta p=eE_{0}/\omega_{0}$ and
$p_{T}=\gamma_{T}m_{e}v_{T}$ is the momentum thermal spread. In this regime,
we find that $\langle 1/\gamma^{3}\rangle$ is not constant, and the electric
field oscillation is instead well described by a triangular shape [22]. We
also performed simulations in the relativistic regime (considering
$E_{0}/(m_{e}c\omega_{0}/e)\simeq 14$ and $\Delta p/m_{e}c=2-800$), where all
conclusions outlined in this Letter for the non relativistic regime hold
qualitatively and quantitatively.

In the relativistic regime, we find also that
$\Gamma_{\textrm{max}}/\omega_{0}\sim 0.1$, where $\omega_{0}\propto\langle
1/\gamma^{3}\rangle^{1/2}$, according to Eq. (2). For waterbag and exponential
distributions of the type $\exp(-\gamma/\gamma_{T})$, we can show that
$\langle 1/\gamma^{3}\rangle=\gamma_{T}^{-1}$ and $(2\gamma_{T})^{-1}$,
respectively, and thus $\omega_{0}\propto\gamma_{T}^{-1/2}$. Hence, the
unstable modes studied in this work are not as severely suppressed for large
$\gamma_{T}$ as streaming instabilities, for which
$\Gamma_{\textrm{max}}\propto\gamma_{T}^{-3/2}$, and may be more easily
excited in extreme astrophysical settings where the Lorentz factor of the
plasma flows is very large. In particular, we find for typical pulsar
parameters (surface field $B\simeq 10^{12}$ G, period $0.1$ s) that
$\omega_{0}\sim 1-10$ GHz [23, 32], such that the instability develops on a
typical time $1/\Gamma\sim 10~{}(p_{T}/\delta p)^{2}/\omega_{0}\sim
10~{}(p_{T}/\delta p)^{2}$ ns for the $p_{T}/\delta p<1$ regime expected in
these scenarios (see Supplemental Material). This should be compared to the
lifetime of the inductive plasma waves, which is well approximated by the time
between pair production bursts, $T_{\textrm{b}}\sim 1$ $\mu$s [33]. We find
that $1/\Gamma<T_{\textrm{b}}$ for $p_{T}/\delta p<10$, which may be achieved
in the plasma trail behind pair production fronts, where the electric field
oscillates inductively, but at an amplitude that is not enough to trigger a
considerable number of pair production events [32]. This instability may also
perturb particle trajectories in inductive waves, previously identified as a
source of linear acceleration emission [34, 35]. We note that the strong
longitudinal magnetic field typical of these environments is not expected to
play a significant role in the development of this instability (contrary to
e.g. filamentation modes), as particle trajectories are purely one-
dimensional.

In conclusion, we have studied the plasma waves parametrically excited in an
inductively oscillatory plasma equilibrium. Our results show that the energy
in the inductive pump wave is transferred to other plasma modes via an
oscillating two-stream instability. Since the pump wave has infinite
wavelength, other fundamental wave depletion mechanisms such as wave breaking
or Landau damping do not operate, regardless of the value of its amplitude
$E_{0}$. We have presented the dispersion relation of these waves for a
waterbag equilibrium distribution function, which captures thermal effects and
is analytically tractable. In general, infinite branches $\omega(k)$ exist,
each branch being a purely real mode when far from other branches in the
$(\omega,k)$ space. However, coupling between different branches yields
unstable modes, with the maximum growth rate
$\Gamma_{\textrm{max}}/\omega_{0}\sim 0.1$. All analytical results have been
confirmed using PIC simulations. We have also investigated the relevance of
this instability in pulsar magnetospheres, and determined that it may be
excited over short distances following pair production bursts in strong
fields. Furthermore, we speculate that it may be relevant in other
astrophysical scenarios, e.g. black hole magnetospheres [36, 37, 38], where
strong rotation-powered electric fields are also self-consistently screened by
plasma currents.

###### Acknowledgements.

F. Cruz and T. Grismayer contributed equally to this work. This work was
supported by the European Research Council (ERC-2015-AdG Grant 695088) and FCT
(Portugal) (grant PD/BD/114307/2016) in the framework of the Advanced Program
in Plasma Science and Engineering (APPLAuSE, FCT grant PD/00505/2012). We
acknowledge PRACE for granting access to MareNostrum, Barcelona Supercomputing
Center (Spain), where the simulations presented in this work were performed.

Supplemental Material

##  I Derivation of the dispersion relation

We start by linearizing the Vlasov equation for positrons and electrons
(superscripts $\pm$, respectively),

$\frac{\partial f_{1}^{\pm}}{\partial t}+v_{\pm}\frac{\partial
f_{1}^{\pm}}{\partial x}\pm\frac{eE_{0}}{m_{e}}\frac{\partial
f_{1}^{\pm}}{\partial v}=\mp\frac{eE_{1}}{m_{e}}\frac{\partial
f_{0}^{\pm}}{\partial v}\ ,$ (S1)

where subscripts $0$ ($1$) correspond to zeroth (first) order, time-varying
quantities. The left hand side of Eq. (S1) can be interpreted as a total time
derivative of $f_{1}^{\pm}$ in phase space. Thus, $f_{1}^{\pm}$ can be
obtained by integrating the right hand side of Eq. (S1),

$f_{1}^{\pm}(x,v,t)=\mp\frac{e}{m_{e}}\int_{-\infty}^{t}\mathrm{d}t^{\prime}E_{1}(x^{\prime},t^{\prime})\frac{\partial
f_{0}^{\pm}(x^{\prime},v^{\prime},t^{\prime})}{\partial v^{\prime}}\ ,$ (S2)

where all primed quantities are taken along the unperturbed trajectories. The
solution in Eq. (S2) can then be used in Poisson’s equation to find the
dispersion relation. Using an anzats of the type
$A_{1}\sim\bar{A}_{1}\exp(i(kx-\omega t))$ for first order quantities, we can
write

$ik\bar{E}_{1}=-\sum_{s=\pm}\frac{4\pi
e^{2}}{m_{e}}\int\mathrm{d}v\int_{-\infty}^{t}\mathrm{d}t^{\prime}\bar{E}_{1}\frac{\partial
f_{0}^{s}}{\partial
v^{\prime}}\exp\left[ik(x^{\prime}_{s}-x_{s})-i\omega(t^{\prime}-t)]\right]\
.$ (S3)

We now attempt to calculate the integrals in Eq. (S3) along the unperturbed
trajectories. These can be written, in the non relativistic regime, as

$v_{\pm}(t)=v_{0}\pm\delta v\cos(\omega_{0}t+\theta)\equiv v_{0}\pm\delta
v\cos\phi\ ,$ (S4)

where the phase $\theta$ is such that $v=v_{0}$ at some reference time
$t_{0}$, and is the same for all particles in the distribution. Eq. (S4)
allows us to write the velocity at any time $t^{\prime}$ as

$v^{\prime}_{\pm}\equiv v_{\pm}(t^{\prime})=v_{0}\pm\delta
v\cos(\omega_{0}(t^{\prime}-t)+\phi)=v_{\pm}\pm\delta
v[\cos(\omega_{0}(t^{\prime}-t)+\phi)-\cos\phi]\ ,$ (S5)

and the difference in positions between times $t^{\prime}$ and $t$ as

$x^{\prime}_{\pm}-x_{\pm}=\int_{t}^{t^{\prime}}\mathrm{d}\tau
v(\tau)=v_{0}(t^{\prime}-t)\pm\frac{\delta
v}{\omega_{0}}[\sin(\omega_{0}(t^{\prime}-t)+\phi)-\sin\phi]\ .$ (S6)

Using the difference in Eq. (S6), we can write the exponential term in Eq.
(S3) for positrons as

$\exp\left[-i(\omega-kv_{0})(t^{\prime}-t)\right]\exp\left[ik\delta
v/\omega_{0}\left(\sin(\omega_{0}(t^{\prime}-t)+\phi)-\sin\phi\right)\right]\
.$ (S7)

For the electron species, this exponential term has the same shape, with
$\delta v\to-\delta v$. This difference does not introduce any change in the
final result, so we proceed with the derivation focusing on the positron
species terms. Let us focus on the last exponential in Eq. (S7). Using the
Bessel identity

$e^{ia\sin x}=\sum_{n=-\infty}^{+\infty}J_{n}(a)e^{inx}\ ,$

where $J_{n}$ is the Bessel function of the first kind and order $n$, we can
write the last exponential in Eq. (S7) as

$\sum_{n=-\infty}^{+\infty}\sum_{m=-\infty}^{+\infty}J_{n}\left(\frac{k\delta
v}{\omega_{0}}\right)J_{m}\left(\frac{k\delta
v}{\omega_{0}}\right)\exp[in\omega_{0}(t^{\prime}-t)]\exp[i(n-m)\phi]\ .$ (S8)

This can be plugged back in Eq. (S3) to obtain

$ik=-\sum_{s}\frac{4\pi
e^{2}}{m_{e}}\int\mathrm{d}v\int_{-\infty}^{t}\mathrm{d}t^{\prime}\frac{\partial
f_{0}^{s}}{\partial
v^{\prime}}\sum_{n,m}J_{n}J_{m}\exp\left[-i(\omega-n\omega_{0}-kv_{0})(t^{\prime}-t)\right]A_{n,m}(\phi)\
,$ (S9)

where we have written $A_{n,m}(\phi)=\exp\left[i(n-m)\phi\right]$ and where we
have omitted the limits of the sums in $n$ and $m$ and the arguments of the
Bessel functions $J_{n}$ and $J_{m}$. We now focus on evaluating the integrals
in Eq. (S9), by considering, for simplicity, the zeroth order distribution
functions as waterbags,

$f_{0}^{\pm}\equiv f_{0}^{\pm}(v,t)=\frac{n_{0}}{\Delta
v}\left[H(v+v_{T}\mp\delta v\cos\phi)-H(v-v_{T}\mp\delta v\cos\phi)\right]\ .$
(S10)

Here, $\Delta v=2v_{T}$ is the thermal spread of the distribution and $H(x)$
is the Heaviside function. The derivative with respect to $v$ of the
distribution function reads $\partial f_{0}^{\pm}/\partial v=n_{0}/\Delta
v(\delta(v+v_{T}\mp\delta v\cos\phi)-\delta(v-v_{T}\mp\delta v\cos\phi))$,
where $\delta(x)=\mathrm{d}H(x)/\mathrm{d}x$.

To evaluate the integral in $v$ in Eq. (S9), we note from Eq. (S4) and (S10),
that the Dirac deltas in $\partial f_{0}/\partial v^{\prime}$ are of the type

$\delta(v^{\prime}+v_{T}\mp\delta
v\cos(\omega_{0}(t^{\prime}-t)+\phi))=\delta(v+v_{T}\mp\delta v\cos\phi)\ ,$
(S11)

and thus $\partial f_{0}/\partial v^{\prime}=\partial f_{0}/\partial v$. We
can finally replace $v_{0}=v\mp\delta v\cos\phi$ for positrons and electrons,
respectively, and the integration in $v$ can be easily performed. The
dispersion relation becomes, at this point,

$\displaystyle ik=$ $\displaystyle-\frac{8\pi e^{2}n_{0}}{m_{e}\Delta
v}\int_{-\infty}^{t}\mathrm{d}t^{\prime}\sum_{n,m}J_{n}J_{m}A_{n,m}(\phi)\times$
$\displaystyle\big{\\{}\exp\left[-i(\omega-n\omega_{0}+kv_{T})(t^{\prime}-t)\right]-\exp\left[-i(\omega-n\omega_{0}-kv_{T})(t^{\prime}-t)\right]\big{\\}}\
,$ (S12)

where we have performed the sum over species. We now change the time
integration variable to $\tau=t-t^{\prime}$. In this case, we have

$\displaystyle ik=$ $\displaystyle\frac{8\pi e^{2}n_{0}}{m_{e}\Delta
v}\int_{0}^{+\infty}\mathrm{d}\tau\sum_{n,m}J_{n}J_{m}A_{n,m}(\phi)\times$
$\displaystyle\big{\\{}\exp\left[i(\omega-n\omega_{0}+kv_{T})\tau\right]-\exp\left[i(\omega-n\omega_{0}-kv_{T})\tau\right]\big{\\}}\
.$ (S13)

The integrals in Eq. (S13) can be readily evaluated, and we find

$1-\sum_{n}J_{n}\frac{\omega_{0}^{2}}{(\omega-n\omega_{0})^{2}-k^{2}v_{T}^{2}}\sum_{m}J_{m}A_{n,m}(\phi)=0\
.$ (S14)

To obtain a dispersion relation that does not depend on $\phi$, we take an
average over one period of the driving oscillation of $A_{n,m}$,

$\langle A_{n,m}(\phi)\rangle=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\phi
A_{n,m}(\phi)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\phi\exp[i(n-m)\phi]=\delta_{n,m}\
,$ (S15)

where we have used the definition of the Kronecker delta $\delta_{n,m}$.
Plugging this in Eq. (S14), the dispersion relation reduces to

$1-\sum_{n=-\infty}^{+\infty}J_{n}^{2}\left(\frac{k\delta
v}{\omega_{0}}\right)\frac{\omega_{0}^{2}}{(\omega-n\omega_{0})^{2}-k^{2}v_{T}^{2}}=0\
.$ (S16)

##  II Mode coupling equations

The dispersion relation obtained analytically and numerically for the
instability described in this manuscript hints that unstable modes are
coupled. Here we derive a closed set of equations in the fluid limit that puts
in evidence this coupling, and use it to interpret why some modes may dominate
in the simulations presented in this work. For simplicity, we present this
derivation in the non relativistic regime. We attribute the subscripts $0$,
$1$ to zeroth order (equilibrium) and first order quantities and the
superscripts $\pm$ to quantities associated with positrons and electrons,
respectively. We start by linearizing the continuity equations,

$\frac{\partial n_{1}^{\pm}}{\partial t}\pm v_{0}(t)\frac{\partial
n_{1}^{\pm}}{\partial x}+n_{0}\frac{\partial v_{1}^{\pm}}{\partial x}=0\ .$
(S17)

Here, $v_{0}(t)$ is the time dependant zeroth order fluid velocity developed
by both species in the oscillating zeroth order electric field $E_{0}(t)$.
Taking the partial time derivative of Eq. (S17), we get

$\frac{\partial^{2}n_{1}^{\pm}}{\partial t^{2}}\pm\frac{\partial
v_{0}}{\partial t}\frac{\partial n_{1}^{\pm}}{\partial x}\pm
v_{0}\frac{\partial^{2}n_{1}^{\pm}}{\partial t\partial
x}+n_{0}\frac{\partial^{2}v_{1}^{\pm}}{\partial t\partial x}=0\ .$ (S18)

Changing the order of the derivatives in the last term, and using the momentum
equation, we have

$\frac{\partial^{2}n_{1}^{\pm}}{\partial t^{2}}\pm\frac{\partial
v_{0}}{\partial t}\frac{\partial n_{1}^{\pm}}{\partial x}\pm
v_{0}\frac{\partial^{2}n_{1}^{\pm}}{\partial t\partial
x}+n_{0}\frac{\partial}{\partial x}\left(-\frac{1}{n_{0}m_{e}}\frac{\partial
p^{\pm}}{\partial x}\pm\frac{e}{m_{e}}E_{1}\mp v_{0}\frac{\partial
v_{1}^{\pm}}{\partial x}\right)=0\ ,$ (S19)

where $p^{\pm}$ is the positron/electron fluid pressure. The last two terms of
Eq. (S19) can be expressed as a function of $n_{1}^{\pm}$ by using Gauss’s law
and the continuity equations, respectively,

$\frac{\partial^{2}n_{1}^{\pm}}{\partial t^{2}}\pm\frac{\partial
v_{0}}{\partial t}\frac{\partial n_{1}^{\pm}}{\partial x}\pm
v_{0}\frac{\partial^{2}n_{1}^{\pm}}{\partial t\partial
x}\pm\omega_{p}^{2}(n_{1}^{+}-n_{1}^{-})-\frac{\gamma^{\pm}T^{\pm}}{m_{e}}\frac{\partial^{2}n_{1}^{\pm}}{\partial
x^{2}}\mp v_{0}\frac{\partial}{\partial x}\left(-\frac{\partial
n_{1}^{\pm}}{\partial t}\mp v_{0}\frac{\partial n_{1}^{\pm}}{\partial
x}\right)=0\ ,$ (S20)

where we have written $\omega_{p}^{2}=4\pi e^{2}n_{0}/m_{e}$ and
$p^{\pm}=\gamma^{\pm}n^{\pm}T^{\pm}$, with $\gamma^{\pm}$ and $T^{\pm}$ being
the adiabatic index and temperature of each fluid, respectively. Eq. (S20) can
be simplified to

$\frac{\partial^{2}n_{1}^{\pm}}{\partial
t^{2}}+\left(v_{0}^{2}-\gamma^{\pm}{v_{\textrm{th}}^{\pm}}^{2}\right)\frac{\partial^{2}n_{1}^{\pm}}{\partial
x^{2}}\pm\frac{\partial v_{0}}{\partial t}\frac{\partial n_{1}^{\pm}}{\partial
x}\pm 2v_{0}\frac{\partial^{2}n_{1}^{\pm}}{\partial t\partial
x}\pm\omega_{p}^{2}(n_{1}^{+}-n_{1}^{-})=0\ ,$ (S21)

with ${v_{\textrm{th}}^{\pm}}^{2}=T^{\pm}/m_{e}$. Eq. (S21) describes two non-
trivially forced and coupled oscillators, $n_{1}^{\pm}$. Assuming now that
$v_{0}(t)=\delta v\cos(\omega_{0}t)=\delta
v/2(\exp(i\omega_{0}t)+\exp(-i\omega_{0}t))$, we can look for wave solutions
to Eq. (S21), which is a generalization of Mathieu’s equation. For a thorough
review of the stability properties of Mathieu’s equation, see ref. [39].

A particularly insightful equation can be obtained by taking the Fourier
transform in space and time of Eq. (S21). After some algebra, we obtain

$A_{0}^{\pm}n_{1}^{\pm}(\omega)\pm A_{1+}n^{\pm}_{1}(\omega+\omega_{0})\pm
A_{1-}n^{\pm}_{1}(\omega-\omega_{0})-A_{2}\left[n^{\pm}_{1}(\omega+2\omega_{0})+n^{\pm}_{1}(\omega-2\omega_{0})\right]=\omega_{p}^{2}n^{\mp}_{1}(\omega)\
,$ (S22)

where $n_{1}^{\pm}(\omega)=n_{1}^{\pm}(\omega,k)$ is the Fourier mode with
frequency $\omega$ and wave vector $k$, and where
$A_{0}^{\pm}=\omega_{p}^{2}+\gamma^{\pm}k^{2}{v_{\textrm{th}}^{\pm}}^{2}-k^{2}\delta
v^{2}/2-\omega^{2}$, $A_{1\pm}=(k\delta v/2)(2\omega\pm\omega_{0})$ and
$A_{2}=(k\delta v/2)^{2}$. Eq. (S22) shows that each mode $\omega$ is coupled
to its neighbours $\omega\pm\omega_{0}$ and $\omega\pm 2\omega_{0}$. This
result is similar to that obtained in other works describing multiple
light/plasma wave interactions originally inspired by ref. [1].

A possible approach to obtain a dispersion relation from Eq. (S22) is to solve
the system of coupled equations for neighbour modes, e.g. $\omega$,
$\omega\pm\omega_{0}$ and $\omega\pm 2\omega_{0}$. This yields a system of
linear equations that can be truncated to any desired neighbour mode order.
This truncation corresponds to an ordering condition in $k\delta
v/\omega_{0}$, since $A_{1\pm}\propto k\delta v/\omega_{0}$ and
$A_{1}\propto(k\delta v/\omega_{0})^{2}$. We have verified numerically that
keeping only modes $\omega$ and $\omega\pm 1$ yields a system of linear
equations whose solution recovers the dispersion relation in Eq. (6), and in
particular the scaling $\Gamma/\omega_{0}\propto(k\delta v/\omega_{0})^{2/3}$.

##  III Scaling of instability growth rate with temperature

We discuss here the scaling with temperature of the growth rate of the
instability presented in this manuscript. For simplicity, we adopt the
normalization $\delta v\to\delta v/c$, $v_{T}\to v_{T}/c$, $k\to
kv_{T}/\omega_{0}$ and $\omega\to\omega/\omega_{0}$. The dispersion relation
in Eq. (4) reads

$1-\sum_{n}J_{n}^{2}\left(k\frac{\delta
v}{v_{T}}\right)\frac{1}{(\omega-n)^{2}-k^{2}}=0\ .$ (S23)

As mentioned in the manuscript, keeping only the terms $n=\pm 1$ in this
series is a reasonable approximation to determine the growth rate of unstable
thermal modes. The solution to Eq. (S23) is, in this approximation, given by

$\omega^{2}=1+F(k)\pm\sqrt{J_{1}^{4}(k\delta v/v_{T})+4F(k)}\ ,$ (S24)

with $F(k)=J_{n}^{2}(k\delta v/v_{T})+k^{2}$. The imaginary component of this
solution, $\Im(\omega)\equiv\Gamma$, is plotted in Figure S1(a) as a function
of $k$ and $\delta v/v_{T}$, and its maximum value for a given $\delta
v/v_{T}$ ratio, $\Gamma_{\textrm{max}}$, is plotted in Figure S1(b) as a
function of $k$. In Figure S1(b) we show also in dashed-dotted and dashed
black lines the analytical estimate in the asymptotic limits of small and
large $\delta v/v_{T}$. These lines are obtained by taking asymptotic limits
of the Bessel function $J_{1}$, with which we can estimate
$\Gamma_{\textrm{max}}\simeq\omega_{0}J_{1}^{2}/2$. For $\delta v/v_{T}\gg 1$,
we have $\Gamma_{\textrm{max}}/\omega_{0}\simeq v_{T}/\pi\delta v$, whereas
for $\delta v/v_{T}\ll 1$ we find
$\Gamma_{\textrm{max}}/\omega_{0}\simeq(\delta v/v_{T})^{2}/8$. Figure S1(b)
also shows the growth rate of the fastest growing modes in several simulations
with different $\delta v/v_{T}$ ratios in black dots, showing a good agreement
with theoretical predictions. All simulations were performed with fixed
$\delta v/c\simeq 0.14$ and varying $v_{T}$. The growth rates were obtained by
plotting the time evolution of the energy in a region in $k$ space in the
vicinity of the fastest growing modes (see Figure 3(a) and (b)), and fitting
exponential functions in the linear phase of the instability.

Figure S1:  Instability growth rate. (a) shows the wavenumber and $\delta
v/v_{T}$ ratio dependence of the growth rate, and (b) illustrates the
dependence of the maximum growth rate for fixed $\delta v/v_{T}$ as a function
of the wavenumber. Dashed-dotted and dashed black lines in (b) represent
theoretical asymptotic limits of the maximum growth rate, whereas black dots
represent simulation results.

## References

  * Nishikawa [1968] K. Nishikawa, J. Phys. Soc. Jpn. 24, 1152 (1968).
  * Drake  _et al._ [1974] J. F. Drake, P. K. Kaw, Y. C. Lee, G. Schmid, C. S. Liu, and M. N. Rosenbluth, Phys. Fluids 17, 778 (1974).
  * Morales  _et al._ [1974] G. J. Morales, Y. C. Lee, and R. B. White, Phys. Rev. Lett. 32, 457 (1974).
  * Freund and Papadopoulos [1980] H. P. Freund and K. Papadopoulos, Phys. Fluids 23, 139 (1980).
  * Tripathi and Liu [1982] V. K. Tripathi and C. S. Liu, Phys. Fluids 25, 629 (1982).
  * Dysthe  _et al._ [1983] K. B. Dysthe, E. Mjølhus, H. L. Pécseli, and K. Rypdal, Phys. Fluids 26, 146 (1983).
  * Kaw and Dawson [1969] P. K. Kaw and J. M. Dawson, Phys. Fluids 12, 2586 (1969).
  * Perkins and Flick [1971] F. W. Perkins and J. Flick, Phys. Fluids 14, 2012 (1971).
  * Porkolab [1974] M. Porkolab, Phys. Fluids 17, 1432 (1974).
  * Bowen  _et al._ [2018] T. A. Bowen, S. Badman, P. Hellinger, and S. D. Bale, Astrophys. J. Lett. 854, L33 (2018).
  * Dorfman and Carter [2016] S. Dorfman and T. A. Carter, Phys. Rev. Lett. 116, 195002 (2016).
  * Matsukiyo and Hada [2003] S. Matsukiyo and T. Hada, Phys. Rev. E 67, 046406 (2003).
  * Edwards  _et al._ [2016] M. R. Edwards, N. J. Fisch, and J. M. Mikhailova, Phys. Rev. Lett. 116, 015004 (2016).
  * Zakharov [1972] V. E. Zakharov, J. Exp. Theor. Phys. 35, 908 (1972).
  * Goldman [1984] M. V. Goldman, Rev. Mod. Phys. 56, 709 (1984).
  * Forslund  _et al._ [1975] D. W. Forslund, J. M. Kindel, and E. L. Lindman, Phys. Fluids 18, 1002 (1975).
  * Weatherall  _et al._ [1981] J. C. Weatherall, M. V. Goldman, and D. R. Nicholson, Astrophys. J. 246, 306 (1981).
  * Thejappa  _et al._ [2012] G. Thejappa, R. J. MacDowall, M. Bergamo, and K. Papadopoulos, Astrophys. J. 747, L1 (2012).
  * Max [1973] C. E. Max, Phys. Fluids 16, 1480 (1973).
  * Note [1] We refer to the oscillating two-stream instability as the instability that develops when two uniform plasma populations (namely electrons and/or positrons) counter-stream with time-dependent average velocities. This is a configuration similar to that in Ref. [1], but in which the pump and all other coupled modes have a similar frequency. This should not be confused with other misleading definitions of the oscillating two-stream instability, discussed in Ref. [15].
  * Qin and Davidson [2014] H. Qin and R. C. Davidson, Phys. Plasmas 21, 064505 (2014).
  * Levinson  _et al._ [2005] A. Levinson, D. Melrose, A. Judge, and Q. Luo, Astrophys. J. 631, 456 (2005).
  * Philippov  _et al._ [2020] A. Philippov, A. Timokhin, and A. Spitkovsky, Phys. Rev. Lett. 124, 245101 (2020).
  * Sturrock [1971] P. A. Sturrock, Astrophys. J. 164, 529 (1971).
  * Ruderman and Sutherland [1975] M. A. Ruderman and P. G. Sutherland, Astrophys. J. 196, 51 (1975).
  * Note [2] The inductive waves and the instability presented in this work also develop for a single oscillating electron population in an ion background. We choose to present the derivation for pair plasmas due to their applicability to pulsar magnetospheres.
  * Bernstein [1958] I. B. Bernstein, Phys. Rev. 109, 10 (1958).
  * Stix [1992] T. H. Stix,  _Waves in Plasmas_ (AIP Press, New York, 1992).
  * Keston  _et al._ [2003] D. A. Keston, E. W. Laing, and D. A. Diver, Phys. Rev. E 67, 036403 (2003).
  * Fonseca  _et al._ [2002] R. A. Fonseca _et al._ , OSIRIS: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators, in  _Computational Science — ICCS 2002_ , edited by P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, and J. J. Dongarra (Springer Berlin Heidelberg, Berlin, Heidelberg, 2002) pp. 342–351.
  * Fonseca  _et al._ [2008] R. A. Fonseca, S. F. Martins, L. O. Silva, J. W. Tonge, F. S. Tsung, and W. B. Mori, Plasma Phys. Control. Fusion 50, 124034 (2008).
  * Cruz  _et al._ [2021] F. Cruz, T. Grismayer, and L. O. Silva, Astrophys. J. 908, 149 (2021).
  * Timokhin [2010] A. N. Timokhin, Mon. Not. R. Astron. Soc. 408, 2092 (2010).
  * Melrose and Luo [2009] D. B. Melrose and Q. Luo, Astrophys. J. 698, 124 (2009).
  * Reville and Kirk [2010] B. Reville and J. G. Kirk, Astrophys. J. 715, 186 (2010).
  * Levinson and Cerutti [2018] A. Levinson and B. Cerutti, Astron. Astrophys. 616, A184 (2018).
  * Chen and Yuan [2020] A. Y. Chen and Y. Yuan, Astrophys. J. 895, 121 (2020).
  * Kisaka  _et al._ [2020] S. Kisaka, A. Levinson, and K. Toma, Astrophys. J. 902, 80 (2020).
  * Kovacic  _et al._ [2018] I. Kovacic, R. Rand, and S. Mohamed Sah, Appl. Mech. Rev. 70, 10.1115/1.4039144 (2018).


# Stability under dwell time constraints:
Discretization revisited

Thomas Mejstrik  and Vladimir Yu. Protasov Submitted: 22. January 2024. This
work was supported by the Austrian Science Foundation (FWF) under Grant
P33352-N. Thomas Mejstrik, University of Vienna<EMAIL_ADDRESS>Vladimir Yu. Protasov, University of L’Aquila, Italy;
<EMAIL_ADDRESS>
###### Abstract

We decide the stability and compute the Lyapunov exponent of continuous-time
linear switching systems with a guaranteed dwell time. The main result asserts
that the discretization method with step size $h$ approximates the Lyapunov
exponent with the precision $C\,h^{2}$, where $C$ is a constant. Let us stress
that without the dwell time assumption, the approximation rate is known to be
linear in $h$. Moreover, for every system, the constant $C$ can be explicitly
evaluated. In turn, the discretized system can be treated by computing the
Markovian joint spectral radius of a certain system on a graph. This gives the
value of the Lyapunov exponent with a high accuracy. The method is efficient
for dimensions up to, approximately, ten; for positive systems, the dimensions
can be much higher, up to several hundreds.

discretization, dynamical system on graphs, extremal norm, joint spectral
radius, linear switching system, stability, Lyapunov exponent, multinorm,
49M25, 93C30, 37C20, 15A60

##  1 Introduction

We consider a linear switching system of the form

$\left\\{\begin{array}[]{l}\dot{\boldsymbol{x}}(t)=A(t){\boldsymbol{x}}(t),\quad
t\in[0,+\infty),\quad A(t)\in{\mathcal{A}}\\\
{\boldsymbol{x}}(0)=0\end{array}\right.$ (1)

with a positive dwell time restriction. This is a linear ODE on the vector-
function ${\boldsymbol{x}}:{\mathbb{R}}_{+}\to{\mathbb{R}}^{d}$ with a matrix
function $A(t)$ taking values from a given finite set
${\mathcal{A}}=\\{A_{1},\ldots,A_{n}\\}$ called _control set_ of matrices
(_regimes, modes_). The _control function_ , or the _switching law_ is an
arbitrary piecewise constant function $A:{\mathbb{R}}_{+}\to{\mathcal{A}}$
with the lengths of every stationary interval (_switching interval_) at least
$m$, where $m>0$ is a given dwell-time constraint. This dwell time assumption
has the practical meaning that the switches between regimes are not
instantaneous but take some positive time.

For the sake of simplicity we consider only the case of finite control sets
${\mathcal{A}}$ and of the same dwell time for all regimes $A_{j}$. All our
results are easily extended to general conditions, see Remark 1.

###  1.1 Statement of the problem

Systems (1) regularly arise in engineering applications, see, for instance,
[2, 16, 17, 32] and references therein. One of the main problems is to find or
estimate the fastest possible growth rate of trajectories, in particular, to
decide about the stability of the system. The stability problem under the
dwell time restrictions have been studied in numerous work [4, 5, 7, 8, 9, 31,
33].

The _L_ yapunov exponent $\sigma=\sigma({\mathcal{A}})$ of a linear switching
system is the infimum of the numbers $\alpha$ such that for every trajectory
${\boldsymbol{x}}(\mathord{\,\cdot\,})$, we have $\|{\boldsymbol{x}}(t)\|\leq
Ce^{\alpha t}$, $t\geq 0$, for some constant $C=C({\boldsymbol{x}})$. Clearly,
if $\sigma<0$, then the system is asymptotically stable, i.e., all its
trajectories tend to zero as $t\to+\infty$. The converse is also true,
although less trivial [4, 21]. Thus, the asymptotic stability is equivalent to
the inequality $\sigma<0$. If, in addition, the system is irreducible, i.e.,
the matrices from ${\mathcal{A}}$ do not share common nontrivial invariant
subspaces, then the inequality $\sigma\leq 0$ is equivalent to the (usual)
stability, when all trajectories are bounded. Let us note that if the matrices
of the control set ${\mathcal{A}}$ have a common invariant subspace, then the
computation of the Lyapunov exponent is reduced to problems in smaller
dimensions by a common matrix factorization [17]. Therefore, in what follows
we additionally assume that the system is irreducible.

Most of results on the stability of a linear switching systems have been
obtained without the dwell time assumption. Usually it is done either by
constructing a Lyapunov function [3, 11, 12, 18, 21], or by approximating of
the trajectories [23, 29, 30]. The latter includes the discretization
approach, which is, of course, the most obvious way to analyse ODEs. Replacing
the derivative $\dot{\boldsymbol{x}}(t)$ by the divided difference
$\frac{{\boldsymbol{x}}(t+h)-{\boldsymbol{x}}(t)}{h}$, we get the Euler
piecewise-linear approximation. Another way of discretization is to replace
the switching law by a piecewise-constant function with intervals being
multiples of a given step size $h>0$. Solving the corresponding ODE in each
interval we obtain ${\boldsymbol{x}}(t+h)=e^{hA_{j}}{\boldsymbol{x}}(t)$,
which gives a piecewise-exponential approximation of the trajectory. Both of
those methods lead to a discrete time switching system of the form
${\boldsymbol{y}}_{k+1}=B(k){\boldsymbol{y}}_{k}$, $k\geq 0$. For the Euler
discretization, the control set ${\mathcal{B}}$ consists of matrices
$B_{j}=I+hA_{j}$, for the second approach, $B_{j}=e^{hA_{j}}$. The Lyapunov
exponent of the discrete system (we denote it by $\sigma_{h}$) can be
efficiently computed by evaluating the _joint spectral radius_ of the matrix
family ${\mathcal{B}}$ [1]. The recent progress in the joint spectral radius
problem [13, 19] allows us to calculate it with a good accuracy or, in most
cases, even to find it precisely. Thus, the Lyapunov exponent
$\sigma({\mathcal{A}})$ of system (1) can be efficiently computed, provided it
is close enough to the Lyapunov exponent $\sigma_{h}$ of its discretization.
The crucial problem is to estimate the precision, i.e., the difference
$|\sigma({\mathcal{A}})-\sigma_{h}|$, depending on the step size $h$.

The discretization method in the stability problem has drawn much attention in
the recent literature and several estimates for the precision
$|\sigma({\mathcal{A}})-\sigma_{h}|$ have been obtained [12, 28, 29, 30]. For
both aforementioned methods, it is linear in $h$, i.e.,
$|\sigma({\mathcal{A}})-\sigma_{h}|\leq Ch$, where, for the constant $C$,
usually only rough upper bounds are known.

###  1.2 Overview of the main results

We establish lower and upper bounds which localize the Lyapunov exponent of
the system (1) to an interval of length at most $C\,h^{2}$, and moreover, show
that $C$ can be explicitly found.

The main result, Theorem 5, states that the second method of discretization,
when $B_{j}=e^{hA_{j}}$, leads to the double inequality
$\sigma_{h}\leq\sigma({\mathcal{A}})\leq\sigma_{h}+C\,h^{2}$, where the
constant $C$ is expressed by means of the so-called _extremal multinorm_. This
multinorm is constructed simultaneously with the computation of $\sigma_{h}$.

The estimate from Theorem 5 gives a very precise method of computation of the
Lyapunov exponent $\sigma({\mathcal{A}})$. The numerical results in dimensions
$d\leq 9$ are given in Section 6. In dimension $d=6$ the precision is usually
between $0.1$ and $0.15$. The computation time on our PC (5 cores, 3.6 GHz) is
about one hour. For positive systems, the method performs much better even in
high dimensions. For $d\leq 200$, the precision does not exceed $0.01$. For
$d\leq 100$, the computation takes a few seconds.

The use of the new estimate, however, is complicated by the fact that the
stability of a discrete system under the dwell time constraint cannot be
analysed by the traditional scheme. For such systems, the Lyapunov function
may not exist at all [10]. The computation of the Lyapunov exponent requires
the concept of _restricted_ or _Markovian_ joint spectral radius [6, 15, 24,
26], which are special cases of the recent theory of dynamical system on
graphs developed in [10, 22, 25]. That is why we need to do preliminary work
to introduce the system on graphs and the concept of Lyapunov multinorm.

###### Remark 1

Our estimates for the approximation rate do not include the number of
matrices, which makes them applicable for arbitrary compact control set
${\mathcal{A}}$. They are also easily generalized to mode-dependent
constraints on the dwell time, when $m$ is a function of the matrix
$A\in{\mathcal{A}}$. Finally, the results can be extended to mixed (discrete-
continuous) systems with hybrid control. In particular, when every switching
from the regime $A_{i}$ to $A_{j}$ is realized by a given linear operator
$E_{ji}$ that can be different from $e^{\,mA_{j}}$.

###  1.3 The structure of the paper

In Section 2 we introduce auxiliary facts and notation such as Markovian joint
spectral radius, dynamical systems, and Lyapunov multinorms. The fundamental
theorem and corollaries are formulated and discussed in Sections 3, the proofs
are given in Section 4. Sections 5 and 6 present the algorithm and analyse
numerical properties.

Throughout the paper we denote vectors by bold letters, $I$ is the identity
matrix, $\rho(X)$ is the spectral radius of the matrix $X$, which is the
maximum modulus of its eigenvalues.

##  2 Preliminary facts. Dynamical systems on graphs

The discretization with step size $h$ approximates the system (1) with the
dwell time parameter $m>0$ by the discrete-time system
${\boldsymbol{y}}_{k+1}=B(k){\boldsymbol{y}}_{k}$, where for each $k\geq 0$,
the matrix $B(k)$ is either $e^{hA_{j}}$ or $e^{mA_{j}}$,
$A_{j}\in{\mathcal{A}}$, see Definition 4 below. The dwell time assumption
imposes a firm restriction to the switching law: It must be a sequence of
blocks of the form $B(k+N)\cdots B(k+1)B(k)=(e^{hA_{j}})^{N}e^{mA_{j}}$, where
$j\in\\{1,\ldots,n\\}$, the length $N\geq 0$ depends on the block, and two
neighbouring blocks must have different modes $j$. Thus, each block has to
begin with $e^{mA_{j}}$ followed by a power of $e^{hA_{j}}$, and then switch
to the next block with a different $j$. The general theory of discrete systems
with constraints on switching laws has been developed in [6, 10, 15, 26] and
then extended to general dynamical systems on graphs. We begin with necessary
definitions and notation.

Consider a directed graph $G$ with vertices $v_{1},\ldots,v_{n}$ and edges
$\ell_{ji}$, which may include loops $\ell_{ii}$. Each vertex $v_{i}$ is
associated to a finite-dimensional linear space $V_{i}$. To every edge
$\ell_{ji}$, if it exists, we associate a linear operator $E_{ji}:V_{i}\to
V_{j}$ and a positive number $h_{ji}$, which is referred to as the _time of
action_ of $E_{ji}$.

###### Definition 2

The _d_ ynamical system on the graph $G$ is the equation
${\boldsymbol{x}}_{k+1}=E(k){\boldsymbol{x}}_{k}$, on the sequence
$\\{{\boldsymbol{x}}_{k}\\}_{k\geq 0}$, where for each $k$, the operator
$E(k)$ is chosen from the set $\\{E_{ji}\\}_{j=1}^{n}$ if
${\boldsymbol{x}}_{k}\in V_{i}$. The point ${\boldsymbol{x}}_{k}$ corresponds
to the time $t_{k}$ and $t_{k+1}=t_{k}+h_{ji}$.

The usual discrete-time linear switching system
${\boldsymbol{x}}_{k+1}=B(k){\boldsymbol{x}}_{k}$,
$B(k)\in{\mathcal{B}}=\\{B_{j}\\}_{j=1}^{n}$, corresponds to the case when $G$
is a complete graph, $V_{j}$ are all equal to ${\mathbb{R}}^{d}$, every vertex
$v_{j}$ has $n$ incoming edges $\ell_{ji},i=1,\ldots,n$, with the operators
$E_{ji}$ equal to the same operator $B_{j}$, and all time intervals $h_{ji}$
are equal to $1$.

Let us have an arbitrary dynamical system on a graph. We assume that each
space $V_{j}$ is equipped with a certain norm $\|\mathord{\,\cdot\,}\|_{j}$.
The collection of those norms $\\{\|\mathord{\,\cdot\,}\|_{j}\\}_{j=1}^{n}$ is
called a _multinorm_. We use the short notation $\|{\boldsymbol{x}}\|$ meaning
that $\|{\boldsymbol{x}}\|=\|{\boldsymbol{x}}\|_{j}$ for ${\boldsymbol{x}}\in
V_{j}$.

Every trajectory $\\{{\boldsymbol{x}}_{k}\\}_{k=0}^{\infty}$ of the system on
$G$ corresponds to an infinite path $v_{i_{0}}\to v_{i_{1}}\to\cdots$, where
${\boldsymbol{x}}_{0}\in V_{i_{0}}$ is a starting point and
${\boldsymbol{x}}_{k+1}=E_{i_{k+1}i_{k}}{\boldsymbol{x}}_{k}$, $k\geq 0$.
Thus, ${\boldsymbol{x}}_{k}\in V_{i_{k}}$ for every $k$. Every point
${\boldsymbol{x}}_{k}$ corresponds to the time
$t_{k}=\sum_{s=1}^{k}h_{i_{s}i_{s-1}}$, which is the total time of the way
from $v_{i_{0}}$ to $v_{i_{k}}$.

The _(Markovian) joint spectral radius_ $\hat{\rho}=\hat{\rho}({\mathcal{A}})$
of the system is

$\hat{\rho}=\lim\limits_{k\to\infty}\max\limits_{\\{{\boldsymbol{x}}_{s}\\}_{s=0}^{k}}\|{\boldsymbol{x}}_{k}\|^{1/t_{k}},$

where the maximum is computed over all trajectories
$\\{{\boldsymbol{x}}_{s}\\}_{s\geq 0}$ with $\|{\boldsymbol{x}}_{0}\|=1$.
Thus, for every trajectory, we have $\|{\boldsymbol{x}}_{k}\|\leq
C\hat{\rho}^{t_{k}}$, $k\in{\mathbb{N}}$. The joint spectral radius is the
rate of the fastest growth of trajectories. See [10] for the correctness of
the definition and for basic properties of this notion.

Let us now consider an arbitrary cycle of the graph $G$: $v_{i_{0}}\to
v_{i_{1}}\to\cdots\to v_{i_{k}}=v_{i_{0}}$. Denote by $\Pi$ the product of
operators $E_{i_{k}i_{k-1}}\cdots E_{i_{1}i_{0}}$ along the cycle. We have
$\Pi{\boldsymbol{x}}_{0}={\boldsymbol{x}}_{k}$. For the spectral radius
$\rho(\Pi)$, which is the maximal modulus of eigenvalues of $\Pi$, we have
${[\rho(\Pi)]^{1/t_{k}}\leq\hat{\rho}}$. Indeed, the left-hand side is the
rate of growth of a periodic trajectory going along that cycle, and it does
not exceed the maximal rate of growth $\hat{\rho}$ over all trajectories. On
the other hand, there are cycles for which the left-hand side is arbitrarily
close to $\hat{\rho}$ [10, 15].

###### Definition 3

A multinorm is called _extremal_ for a system on a graph if for each
$i,j\in\\{1,\ldots,n\\}$ and for every ${\boldsymbol{x}}\in V_{i}$, we have
$\|E_{ji}{\boldsymbol{x}}\|_{j}\leq\hat{\rho}^{h_{ji}}\|{\boldsymbol{x}}\|_{i}$,
provided the edge $\ell_{ji}$ exists.

For the extremal multinorm, for every trajectory, we have
$\|{\boldsymbol{x}}_{k}\|\leq\hat{\rho}^{t_{k}}\|{\boldsymbol{x}}_{0}\|$,
$k\in{\mathbb{N}}$.

The extremal multinorm always exists, provided the system is irreducible. The
reducibility means the existence of subspaces $V_{j}^{\prime}\subset V_{j}$,
$i=1,\ldots,n$, where at least one subspace is nontrivial and at least one
inclusion is strict, such that $E_{ji}V_{i}^{\prime}\subset V_{j}^{\prime}$,
whenever the edge $\ell_{ji}$ exists.

Theorem A [10] _An irreducible dynamical system on a graph possesses an
extremal multinorm_.

Theorem A implies, in particular, that for every irreducible system, we have
$\hat{\rho}>0$. Otherwise, $E_{ji}{\boldsymbol{x}}=0$ for all
${\boldsymbol{x}}\in V_{i}$, i.e., all $E_{ji}$ are equal to zero, in which
case the system is clearly reducible. Therefore, one can always normalize the
operators as $\tilde{E}_{ji}=\hat{\rho}^{-h_{ji}}E_{ji}$, after which the new
dynamical system (on the same graph and with the same time intervals) has the
joint spectral radius equal to one. Moreover, it possesses the same extremal
norm, for which
$\|\tilde{E}_{ji}{\boldsymbol{x}}\|_{j}\leq\|{\boldsymbol{x}}\|_{i}$, hence,
for every trajectory $\tilde{\boldsymbol{x}}_{k}$, the sequence
$\|\tilde{\boldsymbol{x}}_{k}\|$ is non-increasing in $k$.

Now we are going to define the discretization of the dwell time constrained
system (1) and present it as a system on a suitable graph. Then we apply
Theorem A and use an extremal multinorm to estimate the approximation rate.

###### Definition 4

The _$h$ -discretization of the system ${\mathcal{A}}$ with step size $h$_ is
a system ${\mathcal{A}}_{h}$ on a complete graph with $n$ vertices, in which
$V_{j}={\mathbb{R}}^{d}$, $E_{jj}=e^{\,hA_{j}}$, $h_{jj}=h$ for all
$j=1,\ldots,n$, and $E_{ji}=e^{\,mA_{j}}$, $h_{ji}=m$ for all pairs $(i,j)$,
$i\neq j$.

Thus, the $h$-discretization ${\mathcal{A}}_{h}$ is a dynamical system on a
complete graph, where every vertex $v_{j}$ has $n-1$ incoming edges from all
other vertices, all associated to the operator $e^{\,mA_{j}}$ with the time of
action $m$, and also has a loop associated to $e^{\,hA_{j}}$ with time of
action $h$. Respectively, a multinorm
$\|\mathord{\,\cdot\,}\|=\\{\|\mathord{\,\cdot\,}\|_{j}\\}_{j=1}^{n}$ can now
be interpreted as a collection of norms in ${\mathbb{R}}^{d}$, each associated
to the corresponding regime $A_{j}$. See Figure 1 for an example of a system
${\mathcal{A}}$ with three matrices.

The $h$-discretization can also be presented as a discrete-time linear
switching system in ${\mathbb{R}}^{d}$:
${\boldsymbol{x}}_{k+1}=B(k){\boldsymbol{x}}_{k},\allowbreak k\geq 0$, where
the sequence $B(k)$ has values from the control set
$\\{e^{\,mA_{j}},\allowbreak e^{\,hA_{j}}\\}_{j=1}^{n}$ with the following
restriction:

Every element $B(k)$ equal to $e^{\,mA_{j}}$ or $e^{\,hA_{j}}$ can be followed
by either $e^{hA_{j}}$ or $e^{mA_{s}}$, $s\neq j$.

The operators $e^{\,mA_{j}},e^{\,hA_{j}}$ act during the time intervals of
lengths $m$ and $h$ respectively. We denote by $\hat{\rho}({\mathcal{A}}_{h})$
the (Markovian) joint spectral radius of the system ${\mathcal{A}}_{h}$. By
Theorem A, if the family ${\mathcal{A}}$ is irreducible, then for every $h$
its $h$-discretization possesses an extremal norm. Note that this norm can be
different for different $h$.

$V_{1}$$V_{2}$$V_{3}$$e^{mA_{1}}$$e^{mA_{2}}$$e^{mA_{1}}$$e^{mA_{3}}$$e^{mA_{2}}$$e^{mA_{3}}$$e^{hA_{1}}$$e^{hA_{2}}$$e^{hA_{3}}$
Figure 1: The dynamical system on graph, $n=3$.

##  3 The fundamental theorem

Now we are formulating the main result. We consider a linear switching system
${\mathcal{A}}$ given by (1) and its $h$-discretization ${\mathcal{A}}_{h}$.
The joint spectral radius of ${\mathcal{A}}_{h}$ is denoted by
$\hat{\rho}({\mathcal{A}}_{h})$.

###### Theorem 5

Let ${\mathcal{A}}$ be an irreducible continuous-time linear switching system
with the dwell time constraint $m$. Then for every discretization step $h>0$,
we have

$\sigma_{h}\leq\sigma({\mathcal{A}})\leq\sigma_{h}-\frac{1}{m}\ln\left(1-\frac{\|({\mathcal{A}}-\sigma_{h}I)^{2}\|}{8}h^{2}\right),$
(2)

where $\sigma_{h}=\ln\hat{\rho}({\mathcal{A}}_{h})$,
$\|\mathord{\,\cdot\,}\|=\\{\|\mathord{\,\cdot\,}\|_{j}\\}_{j=1}^{n}$ is an
extremal multinorm of ${\mathcal{A}}_{h}$, and

$\|({\mathcal{A}}-\sigma_{h}I)^{2}\|=\max_{j=1,\ldots,n}\|(A_{j}-\sigma_{h}I)^{2}\|_{j}.$

###### Remark 6

The estimate (2) uses two values: the joint spectral radius of the
$h$-discretization $\hat{\rho}({\mathcal{A}}_{h})$ and the operator norms of
$(A_{j}-\sigma_{h}I)^{2}$ in the $j$th component of the extremal multinorm
$\\{\|\mathord{\,\cdot\,}\|_{j}\\}$ for ${\mathcal{A}}_{h}$. They are both
found by the invariant polytope algorithm [10]. We give a brief description in
Section 5.

Writing the Taylor expansion of the logarithm up to the second order, we
obtain the following

###### Corollary 7

Under the assumptions of Theorem 5, we have

$\sigma_{h}\leq\sigma({\mathcal{A}})\leq\sigma_{h}+\frac{\|({\mathcal{A}}-\sigma_{h}I)^{2}\|}{8m}h^{2}+\mathcal{O}(h^{4}),\quad\text{as
$h\to 0$}.$ (3)

###### Remark 8

If $h=m$, then we have basically the discretization of an unrestricted system,
and (3) becomes

$\sigma_{h}\leq\sigma({\mathcal{A}})\leq\sigma_{h}+\frac{\|({\mathcal{A}}-\sigma_{h}I)^{2}\|}{8}h+\mathcal{O}(h^{3}),\quad\text{as
$h\to 0$},$

which again reveals the linear dependence on the discretization step $h$ for
unrestricted systems.

###### Remark 9

If $\left|\left|\mathord{\,\cdot\,}\right|\right|$ is an approximation of an
extremal multinorm of ${\mathcal{A}}_{h}$ up to a factor of $1+\varepsilon$
i.e. (in the notation of Definition 3)

$\hat{\rho}^{h_{ji}}\|{\boldsymbol{x}}\|_{i}\leq\|E_{ji}{\boldsymbol{x}}\|_{j}\leq(1+\varepsilon)\hat{\rho}^{h_{ji}}\|{\boldsymbol{x}}\|_{i},$

then (3) becomes

$\sigma_{h}^{-}\leq\sigma({\mathcal{A}})\leq\sigma_{h}^{+}-\frac{1}{m}\ln\left(1-\frac{\|({\mathcal{A}}-\sigma_{h}^{-}I)^{2}\|}{8}h^{2}\right),$

where $\sigma_{h}^{-}=\ln\hat{\rho}({\mathcal{A}}_{h})$, and
$\sigma_{h}^{+}=\ln(1+\varepsilon)\hat{\rho}({\mathcal{A}}_{h})$.

The quadratic rate of approximation in formula (3) is quite unexpected since
the trajectories of the continuous-time system are approximated by the
trajectories of its $h$-discretization only with the linear rate.
Nevertheless, the approximation of the Lyapunov exponent is quadratic.

###### Remark 10

The performance of the estimate (2) can be spoiled in two cases: either the
dwell time $m$ is too small, or the operator norm of $(A_{j}-\sigma_{h}I)^{2}$
is too large. Note also that (2) is an a posteriori estimate since the
operator norm depends on $h$ and is not known in advance.

Theorem 5 yields the following stability conditions for system (1):

###### Corollary 11

If an $h$-discretization ${\mathcal{A}}_{h}$ is unstable, then so is the
system ${\mathcal{A}}$. If ${\mathcal{A}}_{h}$ is stable and
$\hat{\rho}\leq\bigl{(}1-\frac{h^{2}\|({\mathcal{A}}-(\ln\rho)I)^{2}\|}{8}\bigr{)}^{1/m}$,
then ${\mathcal{A}}$ is stable.

######  Proof 3.12.

The first inequality in (2) implies that if $\hat{\rho}>1$ and hence
$\sigma_{h}>0$, then $\sigma({\mathcal{A}})>0$. The converse is established
similarly.

##  4 Proofs of the main results

The proof of Theorem 5 is based on a simple geometrical argument. If a curve
connects two ends of a segment of length $h$, then the distance from this
curve to the segment does not exceed $h^{2}$ multiplied by the curvature and
by a certain constant. The nontrivial moment is that we need this property in
an arbitrary norm in ${\mathbb{R}}^{d}$ and need to evaluate the constant
depending on this norm. Then we apply this fact to each component
$\|\mathord{\,\cdot\,}\|_{j}$ of the extremal norm of the system
${\mathcal{A}}_{h}$ and estimate the growth of trajectory of the system
${\mathcal{A}}$.

######  Lemma 4.13.

Let $\|\mathord{\,\cdot\,}\|$ be an arbitrary norm in ${\mathbb{R}}^{d}$ and
${\boldsymbol{x}}:[0,h]\to{\mathbb{R}}^{d}$ be a $C^{2}$-curve. Then, for
every $\tau\in[0,h]$, the distance from the point ${\boldsymbol{x}}(\tau)$ to
the segment $[{\boldsymbol{x}}(0),{\boldsymbol{x}}(h)]$ does not exceed
$\frac{h^{2}}{8}\|\ddot{\boldsymbol{x}}\|_{C[0,h]}$.

######  Proof 4.14.

Denote by $\|\cdot\|^{*}$ the dual norm in ${\mathbb{R}}^{d}$, thus
$\|{\boldsymbol{y}}\|^{*}=\sup_{\|{\boldsymbol{x}}\|=1}({\boldsymbol{y}},{\boldsymbol{x}})$.
Let ${\boldsymbol{x}}(0)={\boldsymbol{0}}$ and
${\boldsymbol{x}}(h)={\boldsymbol{a}}$. Let ${\boldsymbol{x}}(\xi)$ be the
most distant point of the arc $\\{{\boldsymbol{x}}(\tau),\tau\in[0,h]\\}$ to
the segment $[{\boldsymbol{0}},{\boldsymbol{a}}]$ and let this maximal
distance be equal to $r$. Suppose ${\boldsymbol{y}}$ is the closest to
${\boldsymbol{x}}(\xi)$ point of that segment and denote
${\boldsymbol{s}}={\boldsymbol{x}}(\xi)-{\boldsymbol{y}}$. Thus,
$\|{\boldsymbol{s}}\|=r$. The segment $[{\boldsymbol{0}},{\boldsymbol{a}}]$
does not intersect the interior of the ball of radius $r$ centred at
${\boldsymbol{x}}(\xi)$. Therefore, by the convex separation theorem, there
exists a linear functional ${\boldsymbol{p}}\in{\mathbb{R}}^{d}$,
$\|{\boldsymbol{p}}\|^{*}=1$, which is non-positive on that segment, non-
negative on the ball, and such that
$({\boldsymbol{p}},{\boldsymbol{s}})=\|{\boldsymbol{s}}\|$. Since the point
${\boldsymbol{y}}$ belongs to the ball, it follows that
$({\boldsymbol{p}},{\boldsymbol{y}})=0$. We have
$({\boldsymbol{p}},{\boldsymbol{s}})=({\boldsymbol{p}},{\boldsymbol{x}}(\xi))-({\boldsymbol{p}},{\boldsymbol{y}})=({\boldsymbol{p}},{\boldsymbol{x}}(\xi))$,
therefore, $({\boldsymbol{p}},{\boldsymbol{x}}(\xi))=r$, see Figure 2.

Defining the function
$f(t)=\bigl{(}{\boldsymbol{p}},{\boldsymbol{x}}(t)\bigr{)}$ we obtain

$f(\xi)-f(0)=\bigl{(}{\boldsymbol{p}},{\boldsymbol{x}}(\xi)-{\boldsymbol{x}}(0)\bigr{)}=\bigl{(}{\boldsymbol{p}},{\boldsymbol{x}}(\xi)\bigr{)}=r.$
(4)

Without loss of generality we assume that $\xi\leq\frac{1}{2}h$, otherwise one
can interchange the ends of the segment $[0,h]$. The Taylor expansion of $f$
at the point $\xi$ gives
$f(t)=f(\xi)+f^{\prime}(\xi)(t-\xi)+\frac{1}{2}f^{\prime}(\eta)(t-\xi)^{2}$,
where $\eta\in[t,\xi]$. The maximum of $f(t)$ is attained at $t=\xi$, hence,
$f^{\prime}(\xi)=0$. For $t=0$, this yields
$f(0)=f(\xi)+\frac{1}{2}f^{\prime\prime}(\eta)\xi^{2}$. Combining with (4), we
obtain

$r=-\frac{1}{2}f^{\prime\prime}(\eta)\xi^{2}=\frac{1}{2}\bigl{(}{\boldsymbol{p}},-\ddot{\boldsymbol{x}}(\eta)\bigr{)}\xi^{2}\leq\frac{1}{2}\|\ddot{\boldsymbol{x}}\|\xi^{2}\leq\frac{h^{2}}{8}\|\ddot{\boldsymbol{x}}\|.$

which completes the proof.

$s$${\boldsymbol{x}}(0)={\boldsymbol{0}}$${\boldsymbol{y}}$${\boldsymbol{x}}(\xi)$${\boldsymbol{x}}(h)={\boldsymbol{a}}$${\boldsymbol{p}}\leq
0$${\boldsymbol{p}}\geq 0$$r\leq\frac{h^{2}}{8}\|\ddot{\boldsymbol{x}}\|$
Figure 2: Construction in proof of Lemma 4.13.

######  Theorem 4.15.

Let ${\boldsymbol{x}}(t)$ be a solution of the differential equation
$\dot{\boldsymbol{x}}=A{\boldsymbol{x}}$ with a constant $d\times d$ matrix
$A$, $\|\mathord{\,\cdot\,}\|$ be an arbitrary norm in ${\mathbb{R}}^{d}$, and
$h\in\bigl{(}0,\sqrt{8/\|A^{2}\|}\bigr{)}$ be a number. Then for every
$\tau\in[0,h]$, we have

$\bigl{\|}{\boldsymbol{x}}(\tau)\bigr{\|}\leq\frac{1}{1-\frac{h^{2}}{8}\|A^{2}\|}\max\Bigl{\\{}\|{\boldsymbol{x}}(0)\|,\|{\boldsymbol{x}}(h)\|\Bigr{\\}}.$
(5)

######  Proof 4.16.

It suffices to consider the vector ${\boldsymbol{x}}(\tau)$ with the maximal
norm over all $\tau\in[0,h]$. Since
$\ddot{\boldsymbol{x}}=A\dot{\boldsymbol{x}}=A^{2}{\boldsymbol{x}}$, it
follows that
$\|\ddot{\boldsymbol{x}}\|_{C[0,h]}=\max_{t\in[0,h]}\|A^{2}{\boldsymbol{x}}(t)\|\leq\|A^{2}\|\cdot\|{\boldsymbol{x}}(\tau)\|$.
Hence, by Lemma 4.13, the distance from the point ${\boldsymbol{x}}(\tau)$ to
the closest point ${\boldsymbol{y}}$ of the segment
$[{\boldsymbol{x}}(0),{\boldsymbol{x}}(h)]$ does not exceed
$\frac{h^{2}}{8}\|A^{2}\|\cdot\|{\boldsymbol{x}}(\tau)\|$. On the other hand,
this distance is not less than
$\|{\boldsymbol{x}}(\tau)\|-\|{\boldsymbol{y}}\|$. It remains to note that
$\|{\boldsymbol{y}}\|\leq\max\\{\|{\boldsymbol{x}}(0)\|,\|{\boldsymbol{x}}(h)\|\\}$,
which follows from the convexity of the norm. Thus,

$\frac{h^{2}}{8}\|A^{2}\|\cdot\|{\boldsymbol{x}}(\tau)\|\geq\|{\boldsymbol{x}}(\tau)\|-\max\Bigl{\\{}\|{\boldsymbol{x}}(0)\|,\|{\boldsymbol{x}}(h)\|\Bigr{\\}}.$

Expressing $\|{\boldsymbol{x}}(\tau)\|$ we arrive at (5).

######  Proof 4.17 (Proof of Theorem 5).

The lower bound follows trivially. To prove the upper bound, we first assume
that $\hat{\rho}({\mathcal{A}}_{h})=1$. By Theorem A, the system
${\mathcal{A}}_{h}$ possesses an extremal multinorm
$\|\mathord{\,\cdot\,}\|=\\{\|\mathord{\,\cdot\,}\|_{j}\\}_{j=1}^{n}$. For
each $j\leq n$, we denote
$\alpha_{j}=-\frac{1}{m}\ln\left(1-\frac{\|A_{j}^{2}\|_{j}}{8}h^{2}\right)$.
Let us show that for every $i\neq j$ and for every point
${\boldsymbol{z}}_{0}\in V_{i}$, the trajectory ${\boldsymbol{z}}(t)$
generated in $V_{j}$ by the ODE $\dot{\boldsymbol{z}}=A_{j}{\boldsymbol{z}}$
and starting at ${\boldsymbol{z}}_{0}$ possesses the property:
$\|{\boldsymbol{z}}(t)\|_{j}\leq e^{\alpha_{j}t}\|{\boldsymbol{z}}_{0}\|_{j}$
for all $t\geq m$. Since the multinorm is extremal, it follows that for every
integer $k\geq 0$, one has
$\bigl{\|}{\boldsymbol{z}}(m+kh)\bigr{\|}_{j}=\bigl{\|}e^{mA_{j}}\bigl{(}e^{hA_{j}}\bigr{)}^{k}{\boldsymbol{z}}_{0}\bigr{\|}_{j}\leq\|{\boldsymbol{z}}_{0}\|_{j}$.
Let $k$ be the maximal integer such that $t\geq m+kh$. Thus, $t=m+kh+\tau$,
$\tau\in[0,h)$. Consider the arc of the trajectory
${\boldsymbol{z}}(\mathord{\,\cdot\,})$ on the time interval $[m+kh,m+(k+1)h]$
and denote ${\boldsymbol{x}}(0)={\boldsymbol{z}}(m+kh)$,
${\boldsymbol{x}}(h)={\boldsymbol{z}}(m+kh+h)$. Observe that both
$\|{\boldsymbol{x}}(0)\|$ and $\|{\boldsymbol{x}}(h)\|$ do not exceed
$\|{\boldsymbol{z}}_{0}\|$. Applying Theorem 4.15 we obtain

$\displaystyle\|{\boldsymbol{z}}(t)\|_{j}$
$\displaystyle=\|{\boldsymbol{x}}(\tau)\|_{j}\leq\frac{1}{1-\frac{h^{2}}{8}\|A_{j}^{2}\|_{j}}\|{\boldsymbol{z}}_{0}\|_{j}$
$\displaystyle=e^{\alpha_{j}m}\|{\boldsymbol{z}}_{0}\|_{j}\leq
e^{\alpha_{j}t}\|{\boldsymbol{z}}_{0}\|_{j}.$

Using the multinorm notation and setting $\alpha=\max_{j=1}^{n}\alpha_{j}$, we
conclude that for every trajectory ${\boldsymbol{z}}(\mathord{\,\cdot\,})$
generated by one regime it holds that $\|{\boldsymbol{z}}(t)\|\leq e^{\alpha
t}\|{\boldsymbol{z}}(0)\|$ for all $t\geq m$.

Consider now an arbitrary trajectory ${\boldsymbol{x}}(\mathord{\,\cdot\,})$
of the system ${\mathcal{A}}$. If it does not have switches, then it is
generated by some regime $A_{j}$, and the proof follows immediately from the
inequality above. Let it have the switching points $t_{0}<t_{1}<\ldots$. This
set can be infinite or finite. Applying the inequality above to arbitrary
$t\in[t_{k},t_{k+1}]$ and denoting
${\boldsymbol{z}}(0)={\boldsymbol{x}}(t_{k})$, we obtain
$\|{\boldsymbol{x}}(t)\|\leq e^{\alpha(t-t_{k})}\|{\boldsymbol{x}}(t_{k})\|$.
Now take a time interval $[0,T]$ and let $t_{N}$ be the largest switching
point on it. Applying our inequality successively for all switching intervals,
we get

$\displaystyle\|{\boldsymbol{x}}(T)\|$ $\displaystyle\leq
e^{\alpha(t-t_{N})}\prod_{k=0}^{N-1}e^{\alpha(t_{k+1}-t_{k})}\|{\boldsymbol{x}}(t_{0})\|$
$\displaystyle\leq e^{\alpha(T-t_{0})}\|{\boldsymbol{x}}(t_{0})\|.$

Thus, for every trajectory ${\boldsymbol{x}}(\mathord{\,\cdot\,})$, we have
$\|{\boldsymbol{x}}(T)\|=\mathcal{O}(e^{\alpha t})$ as $T\to\infty$. Let us
remember that the norm $\|\mathord{\,\cdot\,}\|$ can be different on different
switching intervals. Nevertheless, they all belong to the finite set of norms
$\\{\|\mathord{\,\cdot\,}\|_{j}\\}_{j=1}^{n}$ which are all equivalent.
Therefore, $\sigma({\mathcal{A}})\leq\alpha$.

This concludes the proof for the case $\hat{\rho}=1$. The general case follows
from this one by normalization: We replace the family ${\mathcal{A}}$ by
$\tilde{\mathcal{A}}={\mathcal{A}}-\sigma_{h}I$. Then
$\hat{\rho}(e^{h\tilde{\mathcal{A}}})=1$ and
$\sigma(\tilde{\mathcal{A}})=\sigma({\mathcal{A}})-\sigma_{h}$. Finally,
applying the theorem for the system $\tilde{\mathcal{A}}$ and substituting to
(2), we complete the proof.

Now, to put Theorem 5 into practise, we need to compute the value
$\hat{\rho}({\mathcal{A}}_{h})$ and construct an extremal multinorm. We will
do it in Section 5 for a general system on a graph.

##  5 Computing the joint spectral radius and an extremal multinorm for a
system on a graph

Theorem 5 gives a recipe to compute the Lyapunov exponent of a linear
switching system (1) with sufficiently high accuracy. Due to the quadratic
dependence of $h$ in (3) one can make the distance between the upper and lower
bounds small by choosing an appropriate discretization step $h$. This plan
requires solving two problems: (i) Compute the value of
$\hat{\rho}({\mathcal{A}}_{h})$. (ii) Construct an extremal multinorm for
${\mathcal{A}}_{h}$. Both are solved simultaneously by the invariant polytope
algorithm _(ipa)_ derived in [13]. The ipa computes the joint spectral radius
of several matrices by constructing an extremal norm. In [10] the invariant
polytope algorithm was extended to discrete time systems with restrictions and
to systems of graphs. We present the main idea of the algorithm in this
section; for details see [13, 19, 10, 14]. The reference implementation of the
algorithm can be found at  _gitlab.com/tommsch/ttoolboxes_.

The (Markovian) invariant polytope algorithm (_ipa_) finds the joint spectral
radius $\hat{\rho}({\mathcal{A}}_{h})$ and an extremal multinorm. It consists
of two steps.

#### Step 1.

We fix a number $N$ (not very large) and exhaust all cycles $v_{i_{0}}\to
v_{i_{1}}\to\cdots\to v_{i_{L}}=v_{i_{0}}$ of $G$ of length at most $N$. To
each cycle we denote by $\Pi=E_{i_{L}i_{L-1}}\cdots E_{i_{1}i_{0}}$ the
product of the linear operators associated to its edges and by
$T=\sum_{k=1}^{L}h_{i_{k}i_{k-1}}$ its total time. Recall that
$E_{jj}=e^{hA_{j}}$, $h_{jj}=h$ for all $j$ and $E_{ji}=e^{mA_{j}}$ ,
$h_{ji}=m$ for all pairs $(i,j)$, $i\neq j$.

We choose a cycle with the biggest value $r=\rho(\Pi)^{1/T}$ and call it
_leading cycle_ and respectively  _leading product_ $\Pi$. We set
$\tilde{A}_{j}=A_{j}-(\ln r)I$, $j=1,\ldots,n$. For the new system
$\tilde{\mathcal{A}}_{h}$, we have
$\rho(\tilde{\Pi})=\rho(\tilde{E}_{i_{L}i_{L-1}}\cdots\tilde{E}_{i_{1}i_{0}})=1$.

In many cases to find the leading cycle we can avoid the exhaustion and use
the auxiliary Algorithm .22. It is an adaptation of the modified Gripenberg
algorithm from [19] to our setting. Some of its numerical properties are
assessed in Appendix A.

For the sake of simplicity of exposition, we assume that the leading
eigenvector of $\tilde{\Pi}$ is positive and thus equal to one. The general
case is considered similarly, see [10, 20].

Denoting by ${\boldsymbol{x}}_{0}$ the leading eigenvector of $\tilde{\Pi}$,
it follows that
${\boldsymbol{x}}_{0}\to{\boldsymbol{x}}_{1}\to\cdots\to{\boldsymbol{x}}_{L-1}\to{\boldsymbol{x}}_{0}$
is the periodic trajectory corresponding to that cycle.

#### Step 2.

We try to prove that actually $\hat{\rho}({\mathcal{A}}_{h})=r$ and, if so, to
find an extremal norm.

For every $j$, we denote by ${\mathcal{V}}_{j}^{(0)}$ the set of points
${\boldsymbol{x}}(s)$, $s=0,\ldots,L-1$ that belong to $V_{j}$. If there are
no such points, then ${\mathcal{V}}_{j}^{(0)}=\emptyset$. Suppose after $k$
iterations we have finite sets ${\mathcal{V}}_{j}^{(k)}\subset\allowbreak
V_{j}$, $j=1,\ldots,n$. Denote by
$\operatorname{co}_{s}{{\mathcal{V}}_{j}^{(j)}}$ the convex hull of the set
${\mathcal{V}}_{j}^{(k)}\cup(-{\mathcal{V}}_{j}^{(k)})$. Now for every $j$, we
add to ${\mathcal{V}}_{j}^{(k)}$ all points of the sets
$e^{m\tilde{A}_{j}}{\mathcal{V}}_{s}^{(k)}$ for all $s\neq j$ and of the set
$e^{\,h\tilde{A}_{j}}{\mathcal{V}}_{j}^{(k)}$. We add only those points that
do not belong to $\operatorname{co}_{s}{{\mathcal{V}}_{j}^{(j)}}$, the others
are redundant and we discard them. This way we obtain the sets
${\mathcal{V}}_{j}^{(k+1)}$, $j=1,\ldots,n$. We do this until
${\mathcal{V}}_{j}^{(k+1)}={\mathcal{V}}_{j}^{(k)}$ for all $j$, in which case
the algorithm terminates. We conclude that
$\hat{\rho}({\mathcal{A}}_{h})=\rho(\Pi)^{1/T}$ and the Minkowski norms of the
polytopes $\\{\operatorname{co}_{s}{{\mathcal{V}}_{j}^{(j)}}\\}_{j=1}^{n}$
form an extremal multinorm for ${\mathcal{A}}_{h}$.

######  Remark 5.18.

To find the Lyapunov exponent $\sigma({\mathcal{A}})$, we need to compute the
operator norms $\|A_{j}-\sigma_{h}I\|_{j}$. Since the unit balls of the norms
$\left|\left|\mathord{\,\cdot\,}\right|\right|_{j}$ are polytopes, this can be
done efficiently by solving an LP problem [13].

To achieve a good precision one needs to choose an appropriate step size $h$
which, however, cannot be too small. Otherwise, the matrices $e^{hA_{j}}$ will
be close to the identity complicating the computation of the joint spectral
radius [12]. In particular, the length of the leading cycle may get too large
to be handled efficiently or cannot be found at all, see Example 6.21. In most
cases $h$ cannot be chosen less than $0.1$ (as from our numerical tests).

###  5.1 Positive systems

A continuous-time linear switching system is called _positive_ if every
trajectory ${\boldsymbol{x}}(t)$ starting in the positive orthant
${\mathbb{R}}^{d}_{+}$ remains in ${\mathbb{R}}^{d}_{+}$ for all $t$. Positive
systems have been studied widely in literature due to many applications.

The positivity of the system is equivalent to that all the matrices $A_{j}$
are Metzler, i.e. all off-diagonal entries are nonnegative. In this case all
the matrices $e^{\,hA_{j}},e^{\,mA_{j}}$ of the discrete-time system
${\mathcal{A}}_{h}$ are nonnegative. Hence, the matrix $\Pi$ in the invariant
polytope algorithm is also nonnegative, and therefore, the Perron-Frobenius
theorem implies the nonnegativity of the leading eigenvector
${\boldsymbol{x}}_{0}$. Consequently, all the sets ${\mathcal{V}}^{(k)}_{j}$
are nonnegative, i.e., the algorithm runs entirely in ${\mathbb{R}}^{d}_{+}$.
It is then possible to replace the polytopes
$\operatorname{co}_{s}{{\mathcal{V}}_{j}^{(j)}}$ by the positive polytopes
$\operatorname{co}_{+}{{\mathcal{V}}_{j}^{(j)}}=\\{{\boldsymbol{x}}\geq
0:{\boldsymbol{x}}\leq\operatorname{co}_{s}{{\mathcal{V}}_{j}^{(j)}}\\}$,
where the inequalities are understood element wise.

Since the positive polytopes $Q_{j}^{(k)}$ are in general much larger than
$P_{j}^{(k)}$, they absorb more points in each iteration. This reduces
significantly the number of vertices of the polytopes and, respectively, the
complexity of each iteration. In fact, this modification of the invariant
polytope algorithm for positive systems works very efficiently even for large
dimensions $d$.

##  6 Numerical results

We demonstrate the performance of estimate (2) from Theorem 5 to the
computation of the Lyapunov exponent. Given pairs of random matrices and
random dwell time $m\in(0,1)$, Table 1 presents the performance of the
Markovian ipa. For each dimension $d$, we conducted about 15 tests and we give
the median of the best achieved accuracies i.e. the maximal lower bound from
Equation (3). As test matrices we used $(lhs)$ matrices with normally
distributed entries, and $(rhs)$ Metzler matrices with random integer entries
in $[-9,9]$. To make the examples more interesting, we omit simple cases when
one matrix dominates others (which often occurs). To this end, the matrices
got normalized such that the 2-norm is equal to $1$.

One can see that the Markovian ipa can compute bounds in reasonable time
(although the needed time depends very strongly on the matrices) up to
dimension $\approx 10$ for general matrices, and up to dimension $\approx 500$
for Metzler matrices.

######  Remark 6.19.

The scripts used to obtain the experimental results can be found at
_gitlab.com/tommsch/ttoolboxes/-/tree/master/demo/dwelltime_.

Table 1: Computation of the Lyapunov exponent for arbitrary systems and for
positive systems

Random matrices
---
dim | $\text{ub}-\text{lb}$ | time
$2$ | $0.016512$ |  $4s$
$3$ | $0.020846$ | $58s$
$4$ | $0.036887$ | $780s$
$5$ | $0.108899$ | $2300s$
$6$ | $0.133852$ | $4000s$
$7$ | $0.234899$ | $3300s$
$8$ | $0.314734$ | $4900s$
$9$ | $0.409434$ | $3500s$

Random Metzler matrices
---
dim | $\text{ub}-\text{lb}$ | time
$2$ | $0.018851$ |  $1s$
$5$ | $0.008385$ | $1s$
$13$ | $0.007760$ | $2s$
$34$ | $0.007225$ | $4s$
$89$ | $0.009330$ | $14s$
$144$ | $0.006306$ | $41s$
$233$ | $0.005352$ | $150s$
$377$ | $0.005544$ | $560s$

###  6.1 Examples

Figure 3: Polytopes from Example 6.20.

We begin with a simple two-dimensional example illustrating the Lyapunov
exponent computation by the estimates of Theorem 5 and the Markovian ipa
(Section 5).

######  Example 6.20.

Given dwell time $m=1$, two matrices

$A_{1}=\frac{1}{\sqrt{2}+2}\begin{bmatrix}0&0\\\ 1&0\end{bmatrix},\quad
A_{2}=\frac{1}{\sqrt{2}+2}\begin{bmatrix}-2&-2\\\ -1&-2\end{bmatrix},$

and discretization step $h=0.2$. The product
$\Pi=\allowbreak(e^{mA_{2}})^{1}\allowbreak(e^{hA_{2}})^{7}\allowbreak(e^{mA_{1}})^{1}\allowbreak(e^{hA_{1}})^{180}\allowbreak(e^{mA_{2}})^{1}\allowbreak(e^{hA_{2}})^{7}\allowbreak(e^{mA_{1}})^{1}\allowbreak(e^{hA_{1}})^{181}\allowbreak(e^{mA_{2}})^{1}\allowbreak(e^{hA_{2}})^{7}\allowbreak(e^{mA_{1}})^{1}\allowbreak(e^{hA_{1}})^{183}$
is a leading cycle. The ipa gives $\hat{\rho}({\mathcal{A}}_{h})=\allowbreak
1.0331...$, or respectively $\sigma_{h}=\allowbreak 0.0325...$, and polytopes
$P_{1},\,\allowbreak P_{2}\subseteq\allowbreak{\mathbb{R}}^{2}$ such that
$e^{(m-\sigma_{h})A_{1}}P_{2}\subseteq\allowbreak P_{1}$,
$e^{(h-\sigma_{h})A_{1}}P_{1}\subseteq\allowbreak P_{1}$,
$e^{(h-\sigma_{h})A_{2}}P_{2}\subseteq\allowbreak P_{2}$, and
$e^{(m-\sigma_{h})A_{2}}P_{1}\subseteq\allowbreak P_{2}$.

The relevant norms compute to $\|(A_{1}-\sigma_{h}I)^{2}\|_{P_{1}}=\allowbreak
0.0083...$, and $\|(A_{2}-\sigma_{h}I)^{2}\|_{P_{2}}=\allowbreak 2.8361...$.
Summing up we obtain the bounds
$0.0325<\allowbreak\sigma({\mathcal{A}})<\allowbreak 0.0469$.

In Figure 3 the polytopes $P_{1}$ (blue-thick-dashed line) and $P_{2}$ (red-
thick-dot-dashed line) are plotted, as well the images under the operators
$e^{(m-\sigma_{h})A_{1}}$, …(thin-blue line and thin-red-dot-dashed line).

######  Example 6.21.

Let

$A_{1}\\!=\\!\begin{bmatrix}-1&-1&\phantom{m}1&-1\\\ \phantom{m}1&-1&-1&-1\\\
\phantom{m}1&\phantom{m}1&-1&-1\\\
\phantom{m}1&-1&\phantom{m}1&-1\end{bmatrix},A_{2}\\!=\\!\begin{bmatrix}-1&-1&-1&-1\\\
\phantom{m}1&-1&\phantom{m}1&\phantom{m}1\\\ -1&\phantom{m}1&-1&-1\\\
\phantom{m}1&-1&\phantom{m}1&\phantom{m}1\end{bmatrix},$

and be the dwell time $m=0.5$.

In Table 2 we report the discretization length _( $h$)_, the lower bound
$\sigma_{h}$ of $\sigma({\mathcal{A}})$ given in (2) _(lb)_ , the upper bound
$\sigma_{h}-\frac{1}{m}\ln\left(1-\frac{\|({\mathcal{A}}-\sigma_{h}I)^{2}\|}{8}h^{2}\right)$
of $\sigma({\mathcal{A}})$ _(ub)_ , and a _leading cycle_ , which is a
periodic switching law of the discretized system that produces the fastest
growth of trajectories. The leading cycle is denoted as
$(a^{1}_{1};\,a^{2}_{1};\ldots;\,a^{1}_{k};\,a^{2}_{k})$, where
$a^{1}_{\mathord{\,\cdot\,}}$ and $a^{2}_{\mathord{\,\cdot\,}}$ are the time
of actions of the modes $A_{1}$ and $A_{2}$ respectively. For example, the
leading cycle $(1.3;\,1.7)$ in the first line of Table 2 (the case $h=0.4$) is
the mode $A_{2}$ acting the time $1.3$ followed by the time $1.7$ of the mode
$A_{1}$. From the table we see that the leading cycle seems to stabilize
around $(1.3;\,1.7)$ as $h$ increases.

The character _“?”_ denotes that no leading cycle could be found. On can see,
that for too small values of $h$, the leading cycle cannot be determined, and
thus, the lower bound becomes meaningless. Nevertheless, small discretization
steps may still give good upper bounds.

Table 2: Results for Example 6.21: the bounds for the Lyapunov exponent and the leading cycles depending on $h$.  $h$ | lb | ub | leading cycle
---|---|---|---
$0.400$ | $0.0762$ | $13.2624$ |  $(1.30;\,1.70)$
$0.330$ | $0.0698$ | $8.6828$ | $(1.16;\,1.82;\,1.49;\,1.49)$
$0.300$ | $0.0751$ | $7.1247$ | $(1.40;\,1.70)$
$0.250$ | $0.0742$ | $4.7571$ | $(1.25;\,1.75)$
$0.200$ | $0.0762$ | $3.0066$ | $(1.30;\,1.70)$
$0.125$ | $0.0762$ | $1.1888$ | $(1.25;\,1.625)$
$0.100$ | $0.0000$ | $0.7298$ | $?$
$0.050$ | $0.0000$ | $0.3659$ | $?$
$0.040$ | $0.0000$ | $0.3796$ | $?$
$0.025$ | $0.0000$ | $0.3125$ | $?$

###  .2 Markovian modified Gripenberg algorithm

The Markovian modified Gripenberg algorithm, is an adaption of the modified
Gripenberg algorithm [19]. The main difference is that in each iteration,
still, all possible products have to be considered, but just admissible cycles
(w.r.t. to the graph) are used to compute the intermediate bounds.

######  Algorithm .22 (Markovian modified Gripenberg Algorithm)

The algorithm searches for a leading cycle $\Pi_{\text{max}}$ of a graph $G$
and linear operators $\mathcal{E}$.

$\displaystyle\rho_{\text{max}}=0,\quad\Pi_{0}=I$ $\displaystyle\mathbf{for}\
k=1,\ldots,K$
$\displaystyle\qquad{\mathcal{E}}_{k}=\Pi_{k-1}\times{\mathcal{E}}\quad\text{\footnotesize//
all possible products in iteration $k$}$
$\displaystyle\qquad\rho_{k}=\max\\{\rho(E)^{1/T(E)}:E\in{\mathcal{E}}_{k}\text{
and $E$ is a cycle}\\}$ // where $T(E)$ is the time of action of $E$
$\displaystyle\qquad\Pi_{k}=\\{E\in{\mathcal{E}}_{k}:\rho(E)^{1/T(E)}=\rho_{k}\\}$
$\displaystyle\qquad\mathbf{if}\ \rho_{k}>\rho_{\text{max}}\ \mathbf{then}$
$\displaystyle\qquad\qquad\Pi_{\text{max}}=\Pi_{k},\quad\rho_{\text{max}}=\rho_{k}$
$\displaystyle\qquad\mathbf{else}$
$\displaystyle\qquad\qquad\Pi_{\text{max}}=\Pi_{k}\cup\Pi_{\text{max}}$
$\displaystyle\qquad\Pi_{k}=\Pi_{k}^{+}\cup\Pi_{k}^{-}=\\{E\in{\mathcal{E}}_{k}:\left|\left|E\right|\right|^{1/T(E)}\geq\rho_{k}^{+}\\}\
\cup$
$\displaystyle\qquad\qquad\qquad\\{E\in{\mathcal{E}}_{k}:\left|\left|M\right|\right|^{1/T(M)}\leq\rho_{k}^{-}\\}$
// where $\rho_{k}^{\pm}\geq\rho_{\text{max}}$ and $\rho_{k}^{\pm}$ are such
that $\left|\Pi_{k}^{\pm}\right|=N/2$ $\displaystyle\mathbf{return}\
\Pi_{\text{max}},\rho_{\text{max}}$

######  Example .23.

Given pairs of matrices of various dimensions, random dwell time $m\in(0,1)$
and random discretization step $h\in(0,m)$. Table 3 presents the performance
of the Markovian modified Gripenberg algorithm _(mg)_ compared to a brute
force method _(bf)_. For each dimension $d$ we conducted 15 tests and we
report in how many cases the Gripenberg like algorithm, or the brute force
algorithm worked better Furthermore we give the average runtime of the
algorithm (Note though that the Gripenberg like algorithm in general found the
final result after approximately $2s$).

As test matrices we used $(a)$ matrices with random normally distributed
entries, and $(b)$ Metzler matrices with random integer entries in $[-9,9]$,
always normalized such that the 2-norm equals 1.

One can see that the Gripenberg like algorithm performs in general faster and
better than a brute force method.

Table 3: Results for Example .23: Assessment of Markovian modified Gripenberg
algorithm

$(a)$ Random Gaussian matrices
---
$dim$ | mg better | bf better | $t_{mg}$ |  $t_{bf}$
$2$ | $33\%$ | $\phantom{0}7\%$ | $2.5s$ |  $13.6s$
$3$ | $20\%$ | $\phantom{0}0\%$ | $2.5s$ | $12.7s$
$4$ | $33\%$ | $\phantom{0}7\%$ | $2.4s$ | $16.4s$
$5$ | $13\%$ | $13\%$ | $2.6s$ | $16.1s$
$7$ | $13\%$ | $\phantom{0}0\%$ | $2.6s$ | $12.9s$
$8$ | $20\%$ | $\phantom{0}0\%$ | $2.6s$ | $15.7s$
$11$ | $20\%$ | $\phantom{0}0\%$ | $2.6s$ | $13.3s$
$13$ | $\phantom{0}0\%$ | $\phantom{0}0\%$ | $2.6s$ | $15.9s$
$17$ | $\phantom{0}0\%$ | $\phantom{0}0\%$ | $2.6s$ | $14.0s$
$21$ | $13\%$ | $\phantom{0}0\%$ | $2.6s$ | $17.2s$

$(b)$ Random Metzler matrices
---
$dim$ | mg better | bf better | $t_{mg}$ |  $t_{bf}$
$2$ | $\phantom{0}0\%$ | $0\%$ | $4.2s$ |  $27.3s$
$3$ | $\phantom{0}0\%$ | $0\%$ | $2.6s$ | $27.5s$
$4$ | $40\%$ | $0\%$ | $2.6s$ | $18.6s$
$5$ | $\phantom{0}0\%$ | $0\%$ | $2.6s$ | $26.9s$
$7$ | $20\%$ | $0\%$ | $2.6s$ | $17.6s$
$8$ | $40\%$ | $0\%$ | $2.5s$ | $17.2s$
$11$ | $20\%$ | $0\%$ | $2.6s$ | $13.5s$
$13$ | $\phantom{0}0\%$ | $0\%$ | $2.6s$ | $13.4s$
$17$ | $40\%$ | $0\%$ | $2.5s$ | $13.7s$
$21$ | $20\%$ | $0\%$ | $2.6s$ | $14.2s$

## Acknowledgment

## References

## References

  * [1] N. Barabanov,  Lyapunov indicators of discrete inclusions i–iii,  Autom. Remote Control, 49 (1988), 152–157.
  * [2] C. Basso,  Switch-mode power supplies spice simulations and practical designs,  McGraw-Hill, Inc., New York, NY, USA, 1 edition, 2008.
  * [3] F. Blanchini and S. Miani,  A new class of universal Lyapunov functions for the control of uncertain linear systems,  IEEE Trans. Automat. Control, 44 (1999) 3, 641–647.
  * [4] F. Blanchini, D. Casagrande and S. Miani,  Modal and transition dwell time computation in switching systems: a set-theoretic approach,  Automatica J. IFAC, 46 (2010) 9, 1477–1482.
  * [5] C. Briat and A. Seuret,  Affine characterizations of minimal and mode-dependent dwell times for uncertain linear switched systems,  IEEE Trans. Automat. Control 58 (2013) 5, 1304–1310.
  * [6] X. Dai,  Robust periodic stability implies uniform exponential stability of Markovian jump linear systems and random linear ordinary differential equations,  J. Franklin Inst., 351 (2014), 2910–2937.
  * [7] G. Chesi and P. Colaneri.  Homogeneous rational Lyapunov functions for performance analysis of switched systems with arbitrary switching and dwell time constraints,  IEEE Trans. Automat. Control, 62 (2017) 10, 5124–5137.
  * [8] G. Chesi, P. Colaneri, J. C. Geromel, R. Middleton, and R. Shorten.  A nonconservative LMI condition for stability of switched systems with guaranteed dwell time,  IEEE Trans. Automat. Control, 57 (2012) 5, 1297–1302.
  * [9] Y. Chitour, N. Guglielmi, V. Yu. Protasov, and M. Sigalotti,  Switching systems with dwell time: computing the maximal Lyapunov exponent,  Nonlinear Anal. Hybrid Syst. 40 (2021), Paper No. 101021.
  * [10] A. Cicone, N. Guglielmi, and V. Y. Protasov,  Linear switched dynamical systems on graphs,  Nonlinear Anal. Hybrid Syst., 29 (2018), 165–186.
  * [11] J.C. Geromel and P. Colaneri,  Stability and stabilization of continuous-time switched linear systems,  SIAM J. Control Optim. 45 (2006) 5, 1915–1930.
  * [12] N. Guglielmi, L. Laglia, and V. Yu. Protasov,  Polytope Lyapunov functions for stable and for stabilizable LSS.  Found. Comput. Math., 17 (2017) 2, 567–623.
  * [13] N. Guglielmi and V. Yu. Protasov,  Exact computation of joint spectral characteristics of linear operators,  Found. Comput. Math., 13 (2013), 37–97.
  * [14] V. Yu. Protasov, and R. Kamalov,  Stability of Continuous Time Linear Systems with Bounded Switching Intervals,  SIAM J. Control Optimiz., 61 (2023) 5, 3051–3075.
  * [15] V. Kozyakin,  The Berger-Wang formula for the Markovian joint spectral radius,  Linear Alg. Appl., 448 (2014), 315–328.
  * [16] T. Kröger,  On-line trajectory generation in robotic systems: basic concepts for instantaneous reactions to unforeseen (sensor) events,  Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.
  * [17] D. Liberzon,  Switching in systems and control,  Systems & Control: Foundations & Applications. Birkhäuser Boston, Inc., Boston, MA, 2003.
  * [18] D. Liberzon and A. S. Morse,  Basic problems in stability and design of switched systems,  IEEE Control Systems Magazine, 19 (1999), 59–70.
  * [19] T. Mejstrik,  Improved invariant polytope algorithm and applications, ACM Trans. Math. Softw., ACM Transactions on Mathematical Software 46 (2020) 3, 1–26.
  * [20] T. Mejstrik and V. Yu. Protasov,  Elliptic polytopes and invariant norms of linear operators,  Calcolo 60 (2023) 56.
  * [21] A. Molchanov and Y. Pyatnitskiy,  Criteria of asymptotic stability of differential and difference inclusions encountered in control theory,  Syst. Control Lett., 13 (1989), 59–64.
  * [22] P. Pepe,  Converse Lyapunov theorems for discrete-time switching systems with given switches digraphs, IEEE Trans. Autom. Control, 64 (2019) 6, 2502–2508.
  * [23] A. Pietrus and V. Veliov,  On the discretization of switched linear systems,  Syst. Cont. Letters, 58 (2009) 6, 395–399.
  * [24] M. Philippe, R. Essick, R. Dullerud, and R. M. Jungers,  Stability of discrete-time switching systems with constrained switching sequences,  Automatica, 72 (2016), 242–250.
  * [25] M. Philippe and R. M. Jungers,  A sufficient condition for the boundedness of matrix products accepted by an automaton,  Proceedings of the 18th International Conference on Hybrid Systems: Computation and Control, ACM New York, NY, USA (2015), 51–57.
  * [26] M. Philippe, G. Millerioux, and R. M. Jungers,  Deciding the boundedness and dead-beat stability of constrained switching systems,  Nonlinear Anal. Hybr. Syst., 23 (2017), 287–299.
  * [27] V. Yu.Protasov,  Generalized Markov-Bernstein inequalities and stability of dynamical systems,  Proc. Steklov Inst., 319 (2022), 237–252.
  * [28] V. Yu. Ptotasov and R. M. Jungers,  Analysing the stability of linear systems via exponential Chebyshev polynomials, IEEE Trans. Autom. Control 61 (2016) 3, 795–798.
  * [29] F. Rossi, P. Colaneri, and R. Shorten,  Padé discretization for linear systems with polyhedral lyapunov functions,  IEEE Trans. Autom. Control, 56 (2011) 11, 2717–2722.
  * [30] R. Shorten, M. Corless, S. Sajja and S. Solmaz,  On Padé approximations, quadratic stability and discretization of switched linear systems,  Syst. & Contr. Letters, 60 (2011) 9, 683–689.
  * [31] M. Souza, A. Fioravanti, and R. Shorten,  Dwell-time control of continuous-time switched linear systems,  In Proceedings of the 54th IEEE Conference on Decision and Control (2015), 4661–4666.
  * [32] A. van der Schaft and H. Schumacher,  An introduction to hybrid dynamical systems,  vol. 251 of Lecture Notes in Control and Information Sciences, Springer-Verlag London, Ltd., London, 2000.
  * [33] W. Xiang,  On equivalence of two stability criteria for continuous-time switched systems with dwell time constraint,  Automatica J. IFAC, 54 (2015), 36–40.