# Modulational instability of dust-ion-acoustic waves in pair-ion plasma
having non-thermal non-extensive electrons

M.K. Islam∗,1, A.A. Noman∗∗,1, J. Akter∗∗∗,1, N.A. Chowdhury†,2, A.
Mannan††,1,3,
T.S. Roy‡,4, M. Salahuddin§,1, and A.A. Mamun§§,1 1 Department of Physics,
Jahangirnagar University, Savar, Dhaka-1342, Bangladesh
2 Plasma Physics Division, Atomic Energy Centre, Dhaka-1000, Bangladesh
3 Institut für Mathematik, Martin Luther Universität Halle-Wittenberg, Halle,
Germany
4 Department of Physics, Bangladesh University of Textiles, Tejgaon Industrial
Area, Dhaka, Bangladesh
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The modulational instability (MI) criteria of dust-ion-acoustic (DIA) waves
(DIAWs) have been investigated in a four-component pair-ion plasma having
inertial pair-ions, inertialess non-thermal non-extensive electrons, and
immobile negatively charged massive dust grains. A nonlinear Schrödinger
equation (NLSE) is derived by using reductive perturbation method. The
nonlinear and dispersive coefficients of the NLSE can predict the
modulationally stable and unstable parametric regimes of DIAWs and associated
first and second order DIA rogue waves (DIARWs). The MI growth rate and the
configuration of the DIARWs are examined, and it is found that the MI growth
rate increases (decreases) with increasing the number density of the
negatively charged dust grains in the presence (absence) of the negative ions.
It is also observed that the amplitude and width of the DIARWs increase
(decrease) with the negative (positive) ion mass. The implications of the
results to laboratory and space plasmas are briefly discussed.

###### keywords:

NLSE , Modulational instability , Rogue waves , Pair-ion plasma.

††journal: “Contribution to Plasma Physics”

##  1 Introduction

The ubiquitous existence of pair-ion (PI) in astrophysical environments such
as upper regions of Titan’s atmosphere [1, 2, 3, 4, 5], cometary comae [6],
the ($H^{+}$, $O_{2}^{-}$) and ($H^{+}$, $H^{-}$) plasmas in the D and
F-regions of Earth’s ionosphere [1, 2, 3, 4, 5], and also in the laboratory
experiments namely, neutral beam sources [7], plasma processing reactors [8],
and Fullerene ($C^{+}$, $C^{-}$) [9] have enormously attracted the plasma
physicists to examine nonlinear electrostatic modes in PI plasma. The
properties of the laboratory PI plasma have also been observed by a number of
authors [9, 10, 11]. The dynamics of the PI plasma system and the
configuration of the associated nonlinear electrostatic structures have been
rigorously changed by the presence of negatively charged massive dust grains.
Misra [4] examined dust-ion-acoustic (DIA) shock waves in PI plasma in the
presence of immobile negatively charged dust grains, and reported that the
height of the shock front seems to decrease with the number density of dust
grains. Mushtaq et al. [5] considered a plasma model having inertialess non-
thermal electrons, inertial positive and negative ions, and immobile
negatively charged massive dust grains to investigate DIA solitary waves, and
found that the amplitude of the solitary waves increases with increasing the
number density and charge state of the negative dust grains.

The highly energetic particles associated with long tail in space plasmas have
been identified by Viking spacecraft [12] and Freja satellite [13]. Cairns et
al. [14] proposed non-thermal velocity distribution for these energetic
particles in space plasmas, and observed electrostatic solitary structures. On
the other hand, Tsallis [15] introduced non-extensive $q$-distribution for
explaining the high energy tail in space plasma. The parameter $q$ in the non-
extensive $q$-distribution describes the deviation of the plasma particles
from the thermally equilibrium state. A number of authors [16, 17, 18] have
considered non-thermal non-extensive distribution for investigating nonlinear
electrostatic waves. Dutta and Sahu [18] examined DIA waves (DIAWs) by
considering inertialess non-thermal non-extensive electrons, inertial PI, and
stationary negatively charged dust grains, and highlighted that the angular
frequency of the DIAWs decreases with non-extensivity of electrons but
increases with non-thermality of electrons.

Rogue waves (RWs) are the results of modulational instability (MI) of the
carrier waves, and are governed by the rational solution of standard nonlinear
Schrödinger equation (NLSE) [19, 20, 21], and have been observed in fiber
optics [22], water waves [23], and plasmas [24], etc. Bailung et al. [24]
experimentally observed the RWs in a multi-component plasma, and found that
the amplitude of the RWs is three times greater in comparison with surrounding
normal waves. Abdelwahed et al. [3] observed the effects of the negative ion
on the RWs in PI plasma and determined that the negative ion number density
increases the nonlinearity as well as the amplitude and width of the RWs in PI
plasma. Bains et al. [25] analyzed the MI of ion-acoustic waves (IAWs) and
identified that the critical wave number ($k_{c}$) increases with $q$. El-
Labany et al. [26] investigated the MI of the IAWs and associated RWs in a
three-component plasma system having inertialess iso-thermal electrons and
inertial PI, and found that the amplitude and width of the RWs increase with
the number density and mass of the negative ion. Javidan and Pakzad [27]
considered a three-component dusty plasma containing inertial ions,
inertialess electrons, and immobile dust grains, and examined the instability
of the IAWs, and demonstrated that the instability growth rate decreases with
the increase in the value of negative dust number density. To the best
knowledge of the authors, no attempt has been made to study the MI of DIAWs
and associated DIA RWs (DIARWs) in a four-component plasma having inertialess
non-thermal non-extensive electrons, inertial positive and negative ions, and
stationary negatively charged massive dust grains. The aim of the present
investigation is, therefore, to derive NLSE and investigate MI criteria of
DIAWs in a four-component plasma, and to observe the effects of various plasma
parameters to the formation of first and second order DIARWs.

This manuscript is organized in the following way: The governing equations are
described in Sec. 2. The standard NLSE is derived in Sec. 3. The MI of the
DIAWs and associated DIARWs are examined in Sec. 4. The conclusion of our
present work is provided in Sec. 5.

##  2 Governing Equations

We consider a four-component PI plasma medium having inertial positive ion
(charge $q_{+}=eZ_{+}$ and mass $m_{+}$), inertial negative ion (charge
$q_{-}=-eZ_{-}$ and mass $m_{-}$), inertialess non-thermal non-extensive
electron (charge $q_{e}=-e$ and mass $m_{e}$), and immobile negatively charged
massive dust grains (charge $q_{d}=-eZ_{d}$ and mass $m_{d}$). Overall, the
charge neutrality condition for our plasma model can be written as
$Z_{+}n_{+0}=n_{e0}+Z_{-}n_{-0}+Z_{d}n_{d0}$, where $Z_{+}$ ($Z_{-}$) is the
charge state of the positive (negative) ion, and $Z_{d}$ is the negatively
charged massive dust grains’ charge state. The propagation of the DIAW is
governed by the following equations:

$\displaystyle\frac{\partial n_{+}}{\partial t}+\frac{\partial}{\partial
x}(n_{+}u_{+})=0,$ (1) $\displaystyle\frac{\partial u_{+}}{\partial
t}+u_{+}\frac{\partial u_{+}}{\partial x}=-\frac{\partial\phi}{\partial x},$
(2) $\displaystyle\frac{\partial n_{-}}{\partial t}+\frac{\partial}{\partial
x}(n_{-}u_{-})=0,$ (3) $\displaystyle\frac{\partial u_{-}}{\partial
t}+u_{-}\frac{\partial u_{-}}{\partial x}=\mu_{1}\frac{\partial\phi}{\partial
x},$ (4) $\displaystyle\frac{\partial^{2}\phi}{\partial
x^{2}}=(1-\mu_{2}-\mu_{3})n_{e}-n_{+}+\mu_{2}n_{-}+\mu_{3},\ $ (5)

where $n_{+}$ $(n_{-})$ is the positive (negative) ion number density
normalized by it’s equilibrium value $n_{+0}$ $(n_{-0})$; $u_{+}$ $(u_{-})$ is
the positive (negative) ion fluid speed normalized by the ion-acoustic wave
speed $C_{+}=(Z_{+}k_{B}T_{e}/m_{+})^{1/2}$ with $T_{e}$ being the non-thermal
non-extensive electron temperature and $k_{B}$ being the Boltzmann constant;
$\phi$ is the electrostatic wave potential normalized by $k_{B}T_{e}/e$; the
time and space variables are normalized by
${\omega^{-1}_{p+}}=(m_{+}/4\pi{Z_{+}}^{2}e^{2}n_{+0})^{1/2}$ and
$\lambda_{D+}=(k_{B}T_{e}/4\pi Z_{+}e^{2}n_{+0})^{1/2}$, respectively. Other
plasma parameters can be written as $\mu_{1}=Z_{-}m_{+}/Z_{+}m_{-}$,
$\mu_{2}=Z_{-}n_{-0}/Z_{+}n_{+0}$, and $\mu_{3}=Z_{d}n_{d0}/Z_{+}n_{+0}$. Now,
the expression for the number density of electrons following non-thermal non-
extensive distribution [16] can be written as

$\displaystyle
n_{e}=\big{[}1+A\phi+B\phi^{2}\big{]}\times\big{[}1+(q-1)\phi\big{]}^{\frac{(q+1)}{2(q-1)}},\
$ (6)

where the parameter $q$ stands for the strength of non-extensive system and
the coefficients $A$ and $B$ are defined by
$A=-16q\alpha/(3-14q+15q^{2}+12\alpha)$ and $B=-A(2q-1)$. Here $\alpha$ is a
parameter determining the number of non-thermal electrons in the model.
Williams et al. [17] discussed the range and the validity of ($q$, $\alpha$)
for solitons. In the limiting case $(q\rightarrow 1$ and $\alpha=0$), the
above distribution reduces to the well-known Maxwell-Boltzmann velocity
distribution. For ($q\rightarrow 1$ and $\alpha\neq 0$), the above
distribution reduces to Cairns’ distribution. Now, by substituting Eq. (6)
into Eq. (5) and expanding up to third order of $\phi$, we get

$\displaystyle\frac{\partial^{2}\phi}{\partial
x^{2}}+n_{+}=1+\mu_{2}+\mu_{2}n_{-}+\gamma_{1}\phi+\gamma_{2}\phi^{2}+\gamma_{3}\phi^{3}+\cdot\cdot\cdot,\
$ (7)

where

$\displaystyle\gamma_{1}=[(1-\mu_{2}-\mu_{3})(2A+q+1)]/2,~{}~{}~{}~{}~{}~{}~{}\gamma_{2}=[(1-\mu_{2}-\mu_{3})\\{8B+(q+1)(4A-q+3)\\}]/8,$
$\displaystyle\gamma_{3}=[(1-\mu_{2}-\mu_{3})\\{24B(q+1)-(q+1)(q-3)(6A-3q+5)\\}]/48.\
$

The terms containing $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$ in Eq. (7)
are due to the contribution of the non-thermal non-extensive electrons.

##  3 Derivation of the NLSE

To study the MI of the DIAWs, we want to derive the NLSE by employing the
reductive perturbation method (RPM) and for that case, first we can write the
stretched co-ordinates in the form [28, 29, 30]

$\displaystyle\xi=\epsilon(x-v_{g}t),$ (8) $\displaystyle\tau=\epsilon^{2}t,\
$ (9)

where $v_{g}$ is the group speed and $\epsilon$ ($0<\epsilon<1$) is a small
parameter measuring the weakness of the dispersion. Then, we can write the
dependent variables as [28, 29, 30]

$\displaystyle
n_{+}=1+\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}n_{+l}^{(m)}(\xi,\tau)\exp[il(kx-\omega
t)],$ (10) $\displaystyle
n_{-}=1+\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}n_{-l}^{(m)}(\xi,\tau)\exp[il(kx-\omega
t)],$ (11) $\displaystyle
u_{+}=\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}u_{+l}^{(m)}(\xi,\tau)\exp[il(kx-\omega
t)],$ (12) $\displaystyle
u_{-}=\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}u_{-l}^{(m)}(\xi,\tau)\exp[il(kx-\omega
t)],$ (13)
$\displaystyle\phi=\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}\phi_{l}^{(m)}(\xi,\tau)\exp[il(kx-\omega
t)],\ $ (14)

where $k$ and $\omega$ are the real variables representing the carrier wave
number and frequency, respectively. The derivative operators can be written as

$\displaystyle\frac{\partial}{\partial t}\rightarrow\frac{\partial}{\partial
t}-\epsilon
v_{g}\frac{\partial}{\partial\xi}+\epsilon^{2}\frac{\partial}{\partial\tau},$
(15) $\displaystyle\frac{\partial}{\partial
x}\rightarrow\frac{\partial}{\partial
x}+\epsilon\frac{\partial}{\partial\xi}.\ $ (16)

Now, by substituting Eqs. (8)-(16) into Eqs. (1)-(4) and (7), and collecting
the terms containing $\epsilon$, the first order (when $m=1$ with $l=1$)
reduced equations can be written as

$\displaystyle
u_{+1}^{(1)}=\frac{k}{\omega}\phi_{1}^{(1)},~{}~{}~{}~{}~{}n_{+1}^{(1)}=\frac{k^{2}}{\omega^{2}}\phi_{1}^{(1)},~{}~{}~{}~{}~{}u_{-1}^{(1)}=-\frac{\mu_{1}k}{\omega}\phi_{1}^{(1)},~{}~{}~{}~{}~{}n_{-1}^{(1)}=-\frac{\mu_{1}k^{2}}{\omega^{2}}\phi_{1}^{(1)},\
$

these relations provide the dispersion relation of DIAWs

$\displaystyle\omega^{2}=\frac{k^{2}\big{(}1+\mu_{1}\mu_{2}\big{)}}{k^{2}+\gamma_{1}}.\
$ (17)

Figure 1: The variation of $\omega$ with $k$ for different values of $q$ when
$\alpha=0.5$ (left panel), and the variation of $\omega$ with $k$ for
different values of $\alpha$ when $q=1.5$ (right panel). Other plasma
parameters are $\mu_{1}=0.5$, $\mu_{2}=0.7$, and $\mu_{3}=0.05$.

We have numerically analyzed Eq. (17) to examine the dispersion properties of
DIAWs for different values of $q$ and $\alpha$. The results are displayed in
Fig. 1, which shows that (a) for small wave number, the angular frequency of
the DIAWs exponentially increases and for large wave number, the dispersion
curves become saturated (both left and right panel); (b) the $\omega$
decreases with $q$ (left panel), and this result is a good agreement with the
result of Dutta and Sahu [18]; and (c) as we increase the non-thermality of
the electron then the $\omega$ also increases (right panel), and this result
also coincides with the work of Dutta and Sahu [18]. The second order
equations (when $m=2$ with $l=1$) are given by

$\displaystyle
u_{+1}^{(2)}=\frac{k}{\omega}\phi_{1}^{(2)}+\frac{i(v_{g}k-\omega)}{\omega^{2}}\frac{\partial\phi_{1}^{(1)}}{\partial\xi},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}n_{+1}^{(2)}=\frac{k^{2}}{\omega^{2}}\phi_{1}^{(2)}+\frac{2ik(v_{g}k-\omega)}{\omega^{3}}\frac{\partial\phi_{1}^{(1)}}{\partial\xi},$
$\displaystyle
u_{-1}^{(2)}=-\frac{k\mu_{1}}{\omega}\phi_{1}^{(2)}-\frac{i\mu_{1}(v_{g}k-\omega)}{\omega^{2}}\frac{\partial\phi_{1}^{(1)}}{\partial\xi},~{}~{}~{}~{}~{}n_{-1}^{(2)}=-\frac{\mu_{1}k^{2}}{\omega^{2}}\phi_{1}^{(2)}-\frac{2ik\mu_{1}(v_{g}k-\omega)}{\omega^{3}}\frac{\partial\phi_{1}^{(1)}}{\partial\xi},\
$

and with the compatibility condition, the group velocity of the DIAWs can be
written as

$\displaystyle
v_{g}=\frac{\omega(1-\omega^{2}+\mu_{1}\mu_{2})}{k(1+\mu_{1}\mu_{2})}.\ $ (18)

The variation of $v_{g}$ with $k$ can be seen from the left panel of Fig. 2,
and it is clear from this figure that the group velocity increases with
increasing the value of $\alpha$. Therefore, one can conclude that the group
velocity of DIAWs increases as the non-thermality of electrons increases. The
coefficients of $\epsilon$ when $m=2$ with $l=2$ provides the second order
harmonic amplitudes which are found to be proportional to
$|\phi_{1}^{(1)}|^{2}$

$\displaystyle
n_{+2}^{(2)}=\mu_{4}|\phi_{1}^{(1)}|^{2},~{}~{}~{}~{}~{}u_{+2}^{(2)}=\mu_{5}|\phi_{1}^{(1)}|^{2},~{}~{}~{}~{}~{}n_{-2}^{(2)}=\mu_{6}|\phi_{1}^{(1)}|^{2},~{}~{}~{}~{}~{}u_{-2}^{(2)}=\mu_{7}|\phi_{1}^{(1)}|^{2},~{}~{}~{}~{}~{}\phi_{2}^{(2)}=\mu_{8}|\phi_{1}^{(1)}|^{2},\
$ (19)

where

$\displaystyle\mu_{4}=\frac{3k^{4}+2k^{2}\omega^{2}\mu_{8}}{2\omega^{4}},~{}~{}~{}~{}~{}~{}~{}\mu_{5}=\frac{k^{3}+2k\omega^{2}\mu_{8}}{2\omega^{3}},~{}~{}~{}~{}~{}~{}~{}\mu_{6}=\frac{3k^{4}\mu_{1}^{2}-2k^{2}\omega^{2}\mu_{1}\mu_{8}}{2\omega^{4}},$
$\displaystyle\mu_{7}=\frac{k^{3}\mu_{1}^{2}-2k\omega^{2}\mu_{1}\mu_{8}}{2\omega^{3}},~{}~{}~{}~{}~{}~{}~{}\mu_{8}=\frac{3k^{4}-2\gamma_{2}\omega^{4}-3\mu_{1}^{2}\mu_{2}k^{4}}{2\omega^{2}(4k^{2}\omega^{2}+\gamma_{1}\omega^{2}-k^{2}-k^{2}\mu_{1}\mu_{2})}.\
$

Now, $m=3$ with $l=0$ and $m=2$ with $l=0$ lead to zeroth harmonic modes as
follows:

$\displaystyle
n_{+0}^{(2)}=\mu_{9}|\phi_{1}^{(1)}|^{2},~{}~{}~{}~{}~{}u_{+0}^{(2)}=\mu_{10}|\phi_{1}^{(1)}|^{2},~{}~{}~{}~{}n_{-0}^{(2)}=\mu_{11}|\phi_{1}^{(1)}|^{2},~{}~{}~{}u_{-0}^{(2)}=\mu_{12}|\phi_{1}^{(1)}|^{2},~{}~{}~{}~{}\phi_{0}^{(2)}=\mu_{13}|\phi_{1}^{(1)}|^{2},\
$ (20)

where

$\displaystyle\mu_{9}=\frac{2k^{3}v_{g}+k^{2}\omega+\mu_{13}\omega^{3}}{v_{g}^{2}\omega^{3}},~{}~{}~{}~{}~{}~{}~{}\mu_{10}=\frac{k^{2}+\mu_{13}\omega^{2}}{v_{g}\omega^{2}},~{}~{}~{}~{}~{}~{}~{}\mu_{11}=\frac{2k^{3}v_{g}\mu_{1}^{2}+k^{2}\mu_{1}^{2}\omega-\mu_{1}\mu_{13}\omega^{3}}{v_{g}^{2}\omega^{3}},$
$\displaystyle\mu_{12}=\frac{k^{2}\mu_{1}^{2}-\mu_{1}\mu_{13}\omega^{2}}{v_{g}\omega^{2}},~{}~{}~{}~{}~{}~{}~{}\mu_{13}=\frac{2\gamma_{2}v_{g}^{2}\omega^{3}+k^{2}\mu_{1}^{2}\mu_{2}(2kv_{g}+\omega)-k^{2}(2kv_{g}+\omega)}{\omega^{3}(1+\mu_{1}\mu_{2}-\gamma_{1}v_{g}^{2})}.\
$

Finally, the third harmonic modes, when $m=3$ and $l=1$, with the help of Eqs.
(17)-(20) give a set of equations which can be reduced to the standard NLSE:

$\displaystyle
i\frac{\partial\Phi}{\partial\tau}+P\frac{\partial^{2}\Phi}{\partial\xi^{2}}+Q|\Phi|^{2}\Phi=0,\
$ (21)

where $\Phi=\phi_{1}^{(1)}$, for simplicity. $P$ and $Q$ are the dispersion
and nonlinear coefficients of the NLSE, respectively,

$\displaystyle P=\frac{3v_{g}(v_{g}k-\omega)}{2\omega k},$ (22) $\displaystyle
Q=\frac{2\gamma_{2}\omega^{3}(\mu_{8}+\mu_{13})-k^{2}\omega(\mu_{4}+\mu_{9}+\mu_{1}\mu_{6}+\mu_{1}\mu_{11})+3\gamma_{3}\omega^{3}-2k^{3}(\mu_{5}+\mu_{10}+\mu_{1}\mu_{7}+\mu_{1}\mu_{12})}{2k^{2}(1+\mu_{1}\mu_{2})}.\
$ (23)

The space and time evolution of the DIAWs in PI plasma are directly governed
by the coefficients $P$ and $Q$, and indirectly governed by the different
plasma parameters such as $q$, $\alpha$, $\mu_{1}$, $\mu_{2}$, and $\mu_{3}$.

Figure 2: The variation of $v_{g}$ with $k$ for different values of $\alpha$
when $q=1.5$ (left panel), and the variation of $P/Q$ with $k$ for different
values of $q$ when $\alpha=0.5$ (right panel). Other plasma parameters are
$\mu_{1}=0.5$, $\mu_{2}=0.7$, and $\mu_{3}=0.05$. Figure 3: The variation of
$\Gamma_{g}$ with $k_{M}$ for different values of $\mu_{3}$ when $q=1$ and
$\alpha=0$ (left panel), and the variation of $\Gamma_{g}$ with $k_{M}$ for
different values of $\mu_{3}$ when $q=1.5$ and $\alpha=0.5$ (right panel).
Other plasma parameters are $k=0.7$, $\phi_{0}=0.5$, $\mu_{1}=0.5$, and
$\mu_{2}=0.7$. Figure 4: The variation of $\Gamma_{g}$ with $k_{M}$ for
different values of $\mu_{3}$ when $\mu_{2}=0$ (left panel), and the variation
of $\Gamma_{g}$ with $k_{M}$ for different values of $\mu_{3}$ when
$\mu_{2}=0.5$ (right panel). Other plasma parameters are $k=0.7$,
$\phi_{0}=0.5$, $\alpha=0.5$, $q=1.5$, and $\mu_{1}=0.5$. Figure 5: The
variation of $|\Phi_{1}|$ with $\xi$ for different values of $\mu_{2}$ when
$q=1$ and $\alpha=0$ (left panel), and the variation of $|\Phi_{1}|$ with
$\xi$ for different values of $\mu_{2}$ when $q=1.5$ and $\alpha=0.5$ (right
panel). Other plasma parameters are $k=0.7$, $\tau=0$, $\mu_{1}=0.5$, and
$\mu_{3}=0.05$. Figure 6: The variation of $|\Phi_{1}|$ with $\xi$ for
different values of $\mu_{1}$ (both left and right panel). Other plasma
parameters are $k=0.7$, $\tau=0$, $\alpha=0.5$, $q=1.5$, $\mu_{2}=0.7$, and
$\mu_{3}=0.05$. Figure 7: The variation of $|\Phi_{2}|$ with $\xi$ for
different values of $\mu_{1}$ (both left and right panel). Other plasma
parameters are $k=0.7$, $\tau=0$, $\alpha=0.5$, $q=1.5$, $\mu_{2}=0.7$, and
$\mu_{3}=0.05$. Figure 8: The variation of first-order (dashed green curve)
and second-order (solid red curve) rational solutions of NLSE at $k=0.7$,
$\tau=0$, $\alpha=0.5$, $q=1.5$, $\mu_{1}=0.5$, $\mu_{2}=0.7$, and
$\mu_{3}=0.05$.

##  4 Modulational Instability and Rogue Waves

The stable and unstable parametric regimes of the DIAWs are determined by the
sign of the dispersion ($P$) and nonlinear ($Q$) coefficients of the standard
NLSE [31, 32, 33, 34]. When $P$ and $Q$ have same sign (i.e., $P/Q>0$), the
evolution of the DIAWs amplitude is modulationally unstable. On the other
hand, when $P$ and $Q$ have opposite sign (i.e., $P/Q<0$), the DIAWs are
modulationally stable in the presence of external perturbations. The plot of
$P/Q$ against $k$ yields stable and unstable parametric regimes of DIAWs. The
point, at which transition of $P/Q$ curve intersects with $k$-axis, is known
as threshold or critical wave number $k$ ($=k_{c}$) [31, 32, 33, 34]. We have
numerically analyzed the variation of $P/Q$ with $k$ for different values of
$q$ in the right panel of Fig. 2, and it can be seen from this figure that (a)
both modulationally stable and unstable parametric regimes of DIAWs can be
obtained; (b) the modulationally stable parametric regime of the DIAWs
increases with an increase in the value of $q$, and this result is a good
agreement with the result of Bains et al. [25]; and (c) the modulationally
unstable parametric regime allows to generate first and second order DIARWs.
When $P/Q>0$ and $k_{M}<k_{c}$, the MI growth rate $(\Gamma_{g})$ is given by
[34, 35]

$\displaystyle\Gamma_{g}=|P|k_{M}^{2}\sqrt{\frac{k_{c}^{2}}{k_{M}^{2}}-1},$
(24)

where $k_{M}$ is the modulated wave number. Now, we have graphically shown how
$\Gamma_{g}$ varies with $k_{M}$ for different values of $\mu_{3}$ in Fig. 3.
It is obvious from Fig. 3 that (a) the growth rate initially increases with
$k_{M}$ and becomes maximum, then reduces to zero; (b) under the consideration
of iso-thermal electrons (i.e., $q=1$ and $\alpha=0$), the maximum value of
the growth rate increases with the number density of the negatively charged
dust grains when other plasma parameters are constant (left panel); (c)
$\Gamma_{g}$ as well as nonlinearity of the plasma system increases
(decreases) with increasing the number density of negative (positive) dust
grains (ion) for their constant charge state under the consideration of non-
thermal non-extensive electrons (right panel); (d) the maximum value of the
growth rate increases rigorously in the presence of non-thermal non-extensive
electrons than the presence of iso-thermal electrons.

The existence of the negative ions in the PI plasma rigourously changes the
dynamics of the system. Fig. 4 describes the variation of $\Gamma_{g}$ with
$k_{M}$ for different values of $\mu_{3}$. It is obvious from this figure that
(a) when we neglect the contribution of negative ions from our plasma system
that means only inertial positive ions, inertialess non-thermal non-extensive
electrons, and immobile dust grains are present in the plasma system, then the
maximum value of $\Gamma_{g}$ decreases with the increase of negative dust
grain number density, and similar effect has also been observed by Javidan and
Pakzad [27] in a three-component plasma system (left panel); (b) but when we
include negative ions in the plasma system, that means a four-component PI
plasma system, then $\Gamma_{g}$ increases with increasing the negatively
charged dust grains (right panel). So, the quasi-neutrality condition and
associated dynamics of the plasma system are fully changed by negative ions.

The first-order rational solution of Eq. (21), which can predict the
concentration of large amount of energy in a small region of the
modulationally unstable parametric regime ($P/Q>0$) of DIAWs, can be written
as [36, 37]

$\displaystyle\Phi_{1}(\xi,\tau)=\sqrt{\frac{2P}{Q}}\bigg{[}\frac{4(1+4iP\tau)}{1+16P^{2}\tau^{2}+4\xi^{2}}-1\bigg{]}\exp(2iP\tau).\
$ (25)

We have plotted Eq. (25) in Figs. 5 and 6 to understand the nonlinear
properties of the PI plasma system as well as the mechanism to the formation
of DIARWs associated with DIAWs in the modulationally unstable parametric
regime. Figure 5 indicates that (a) the amplitude and width of the DIARWs
increase with an increase in the value of $\mu_{2}$ (both left and right
panel); (b) under the consideration of iso-thermal electrons, as we increase
(decrease) the value of negative (positive) ion number density, the amplitude
and width of the DIARWs increase (decrease) when their charge state remain
constant (left panel), and this result is similar with the result of El-Labany
et al. [26]; (c) similarly, the amplitude and width of the DIARWs increase
with increasing (decreasing) negative (positive) ion number density when the
plasma system has non-extensive non-thermal electrons, and this result
coincides with the result of Abdelwahed et al. [3]; (d) the direction of the
variation of amplitude and width has not been changed with $\mu_{2}$ under the
consideration of iso-thermal or non-thermal non-extensive electrons, but the
variation of amplitude and width is severely happened for non-thermal non-
extensive electrons than iso-thermal electrons.

It can be seen from the literature that the PI plasma system can support these
conditions: $m_{-}>m_{+}$ (i.e., $H^{+}$, $O_{2}^{-}$) [2, 3], $m_{-}=m_{+}$
(i.e., $H^{+}$, $H^{-}$) [2, 3], and $m_{-}<m_{+}$ (i.e., $Ar^{+}$, $F^{-}$)
[2, 3]. So, in our present investigation, we have graphically observed the
variation of the electrostatic wave potential with $\mu_{1}$ in Fig. 6 under
the consideration of $m_{-}>m_{+}$ (i.e., $\mu_{1}<1$ and left panel) and
$m_{-}<m_{+}$ (i.e., $\mu_{1}>1$ and right panel), and it is obvious from this
figure that (a) the shape of the first order DIARWs is so much sensitive to
the change of mass and charge state of the positive and negative ions; (b) the
electrostatic wave potential increases with negative ion mass but decreases
with the increase of positive ion mass when their charge state remain constant
(both left and right panel); (c) physically, the nonlinearity as well as the
amplitude of the DIARWs increases (decreases) with negative (positive) ion
mass, and this result is a good agreement with the result of El-Labany et al.
[26].

The second-order rational solution of NLSE can be written as [36, 37]

$\displaystyle\Phi_{2}(\xi,\tau)=\sqrt{\frac{P}{Q}}\bigg{[}1+\frac{G_{2}+iM_{2}}{T_{2}}\bigg{]}\exp(iP\tau),\
$ (26)

where

$\displaystyle
G_{2}(\xi,\tau)=\frac{-\xi^{4}}{2}-6(P\xi\tau)^{2}-10(P\tau)^{4}-\frac{3\xi^{2}}{2}-9(P\tau)^{2}+\frac{3}{8},$
$\displaystyle
M_{2}(\xi,\tau)=-P\tau\Big{[}\xi^{4}+4(P\xi\tau)^{2}+4(P\tau)^{4}-3\xi^{2}+2(P\tau)^{2}-\frac{15}{4}\Big{]},$
$\displaystyle
D_{2}(\xi,\tau)=\frac{\xi^{6}}{12}+\frac{\xi^{4}(P\tau)^{2}}{2}+\xi^{2}(P\tau)^{4}+\frac{\xi^{4}}{8}+\frac{9(P\tau)^{4}}{2}-\frac{3(P\xi\tau)^{2}}{2}+\frac{9\xi^{2}}{16}+\frac{33(P\tau)^{2}}{8}+\frac{3}{32}.\
$

We have graphically shown Eq. (26) in Fig. 7. The nature of the second order
DIARWs is very sensitive to the change of mass and charge state of the PI
under the consideration of $m_{-}>m_{+}$ (i.e., $\mu_{1}<1$ and left panel)
and $m_{-}<m_{+}$ (i.e., $\mu_{1}>1$ and right panel). The amplitude of the
second order DIARWs associated with DIAWs in the modulationally unstable
parametric regime increases (decreases) with the rise of the mass of the
negative (positive) ion for a fixed value of their charge state (both left and
right panel). The physics of this result is that the nonlinearity of the
plasma increases with the mass of negative ion but decreases with the mass of
positive ion.

The comparison of the first and second order DIARWs can be observed from Fig.
8, and it is clear from this figure that (a) the second order DIARWs has
double structures compared with the first-order DIARWs; (b) the amplitude of
the second order DIARWs is always greater than the amplitude of the first
order DIARWs; (c) the potential profile of the second order DIARWs becomes
more spiky (i.e., the taller amplitude and narrower width) than the first
order DIARWs; (d) the second (first) order DIARWs has four (two) zeros
symmetrically located on the $\xi$-axis; (e) the second (first) order DIARWs
has three (one) local maxima.

##  5 Conclusion

We have studied an unmagnetized realistic PI plasma system having inertialess
non-thermal non-extensive electrons, inertial negative and positive ions, and
immobile negatively charged dust grains. The RPM is used to derive the NLSE,
and the nonlinear and dispersive coefficients of the NLSE determine the stable
and unstable parametric regimes of DIAWs. The results that have been found
from our investigation can be summarized as follows:

  * 1.

The angular frequency of the DIAWs decreases (increases) with the increase in
the value of non-extensivity (non-thermality) of the electrons.

  * 2.

Both modulationally stable and unstable parametric regimes of DIAWs have been
observed.

  * 3.

The modulational instability growth rate increases (decreases) with increasing
the number density of the negatively charged dust grains in the presence
(absence) of the negative ions.

  * 4.

The amplitude and width of the DIARWs associated with DIAWs increase
(decrease) with the negative (positive) ion mass.

It may be noted here that the gravitational and magnetic effects are very
important but beyond the scope of our present work. In future and for better
understanding, someone can investigate the nonlinear propagation in a four-
component PI plasma by considering the gravitational and magnetic effects. The
results of our present investigation will be useful in understanding the
nonlinear phenomena in both astrophysical environments such as upper regions
of Titan’s atmosphere [1, 2, 3, 4, 5], cometary comae [6], the ($H^{+}$,
$O_{2}^{-}$) and ($H^{+}$, $H^{-}$) plasmas in the D and F-regions of Earth’s
ionosphere [1, 2, 3, 4, 5], and also in laboratory plasmas namely, neutral
beam sources [7], plasma processing reactors [8], and Fullerene ($C^{+}$,
$C^{-}$) [9], etc.

## Acknowledgements

Authors would like to acknowledge “UGC research project 2018-2019” for their
financial supports to complete this work.

## References

  * [1] H. Massey, Negative Ions, 3rd ed., Cambridge University Press, Cambridge, (1976).
  * [2] R. Sabry, W. M. Moslem, P. K. Shukla, Phys. Plasmas 16, 032302 (2009).
  * [3] H. G. Abdelwahed, E. K. El-Shewy, M. A. Zahran, S. A. Elwakil, Phys. Plasmas 23, 022102 (2016).
  * [4] A. P. Misra, Physics of Plasmas 16, 033702 (2009).
  * [5] A. Mushtaq, M. N. Khattak, Z. Ahmad, and A. Qamar, Phys. Plasmas 19, 042304 (2012).
  * [6] P. H. Chaizy, H. Rème, J. A. Sauvaud, C. D’uston, R. P. Lin, D. E. Larson, et al., Nature (London) 349, 393 (1991).
  * [7] M. Bacal and G. W. Hamilton, Phys. Rev. Lett. 42, 1538 (1979).
  * [8] R. A. Gottscho and C. E. Gaebe, IEEE Trans. Plasma Sci. 14, 92 (1986).
  * [9] W. Oohara and R. Hatakeyama, Phys. rev. lett. 91, 205005 (2003).
  * [10] R. Hatakeyama and W. Oohara, Phys. Scripta 116, 101 (2005).
  * [11] W. Oohara, D. Date, and R. Hatakeyama, Phys. rev. lett. 95, 175003 (2005).
  * [12] R. Boström, IEEE Trans. Plasma Sci. 20, 756 (1992).
  * [13] P. O. Dovner, A. I. Eriksson, R. Boström, and B. Holback, Geophys. Res. Lett. 21, 1827 (1994).
  * [14] R. A. Cairns, A. A. Mamun, R. Bingham, R. Boström, R. O. Dendy, et al., Geophys. Res. Lett. 22, 2709 (1995).
  * [15] C. Tsallis, Journal of statistical physics 52, 479 (1988).
  * [16] M. Tribeche, R. Amour, and P. K. Shukla, Phys. Rev. E 85, 037401 (2012).
  * [17] G. Williams, I. Kourakis, F. Verheest, and M. A. Hellberg, Phys. Rev. E 88, 023103 (2013).
  * [18] D. Dutta and B. Sahu, Commun. Theor. Phys. 68, 117 (2017).
  * [19] N. Akhtar, S. A. El-Tantawy, S. Mahmood, and AM. Wazwaz, Rom. Rep. Phys. 71, 403 (2019).
  * [20] B. Sun and AM. Wazwaz, Comm. Nonlinear Sci. Numer. Simulat. 64, 1 (2018).
  * [21] Shalini, A. P. Misra, and N. S. Saini, J Theor. Appl. Phys. 11, 217 (2017).
  * [22] B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, et al., Nat. Phys. 6, 790 (2010).
  * [23] A. Chabchoub, N. P. Hoffmann, and N. Akhmediev, Phys. Rev. Lett. 106, 204502 (2011).
  * [24] H. Bailung, S. K. Sharma, and Y. Nakamura, Phys. Rev. Lett. 107, 255005 (2011).
  * [25] A. S. Bains, M. tribeche, T. S. Gill, Phys. Plasmas 18, 022108 (2011).
  * [26] S. K. El-Labany, W. M. Moslem, N. A. El-Bedwehy, et al., Astrophys. Space Sci. 338, 3 (2012).
  * [27] K. Javidan, H. R. Pakzad, Astrophys Space Sci 350, 557 (2014).
  * [28] S. Jahan, et al., Plasma Phys. Rep. 46, 90 (2020); N. A. Chowdhury, et al., Chaos 27, 093105 (2017); N. Ahmed, et al., Chaos 28, 123107 (2018); N. A. Chowdhury, et al., Vacuum 147, 31 (2018); M. Hassan, et al., Commun. Theor. Phys. 71, 1017 (2019); R. K. Shikha, et al., Eur. Phys. J. D 73, 177 (2019); N. A. Chowdhury, et al., Plasma Phys. Rep. 45, 459 (2019); T. I. Rajib, et al., Phys. Plasmas 26, 123701 (2019).
  * [29] M. H. Rahman, N. A. Chowdhury, A. Mannan, M. Rahman, and A. A. Mamun, Chinese J. Phys. 56, 2061 (2018).
  * [30] M. H. Rahman, A. Mannan, N. A. Chowdhury, and A. A. Mamun, Phys. Plasmas 25, 102118 (2018).
  * [31] I. Kourakis and P.K. Shukla, Phys. Plasmas 10, 3459 (2003).
  * [32] I. Kourakis and P.K. Shukla, Nonlinear Proc. Geophys. 12, 407 (2005).
  * [33] R. Fedele, Phys. Scr. 65, 502 (2002).
  * [34] S. Sultana and I. Kourakis, Plasma Phys. Control. Fusion 53, 045003 (2011).
  * [35] S. Jahan, et al., Commun. Theor. Phys. 71, 327 (2019); N. A. Chowdhury, et al., Phys. plasmas 24, 113701 (2017); N. A. Chowdhury, et al., Contrib. Plasma Phys. 58, 870 (2018); S.K. Paul, et al., Pramana J. Phys. 94, 58 (2020).
  * [36] A. Ankiewicz, N. Devine, and N. Akhmediev, Phys. Lett. A 373, 3997 (2009).
  * [37] S. Guo, L. Mei and A. Sun, Ann. Phys. 332, 38 (2012).


# Topological Linear System Identification
via Moderate Deviations Theory

Wouter Jongeneel, Tobias Sutter and Daniel Kuhn The authors are with the Risk
Analytics and Optimization Chair, École Polytechnique Fédérale de Lausanne,
{wouter.jongeneel, tobias.sutter<EMAIL_ADDRESS>This research was
supported by the Swiss National Science Foundation under the NCCR Automation,
grant agreement 51NF40_180545.

###### Abstract

Two dynamical systems are topologically equivalent when their phase-portraits
can be morphed into each other by a homeomorphic coordinate transformation on
the state space. The induced equivalence classes capture qualitative
properties such as stability or the oscillatory nature of the state
trajectories, for example. In this paper we develop a method to learn the
topological class of an unknown stable system from a single trajectory of
finitely many state observations. Using tools from moderate deviations theory
we prove that the probability of misclassification decays exponentially with
the number of observations at a rate that is proportional to the square of the
smallest singular value of the true system matrix.

##  I Introduction

We consider the discrete-time linear time-invariant system

$x_{t+1}=\theta x_{t}+w_{t},\quad x_{0}\sim\nu,$ (1)

where $x_{t}\in\mathbb{R}^{n}$ and $w_{t}\in\mathbb{R}^{n}$ denote the state
and the exogenous noise at time $t\in\mathbb{N}$, while $\theta$ represents
the system matrix, and $\nu$ stands for the marginal distribution of the
initial state $x_{0}$. Except for asymptotic stability we assume that nothing
is known about $\theta$, and we aim to identify $\theta$ from a single
trajectory of states $\\{\widehat{x}_{t}\\}_{t=0}^{T}$ generated by (1). A
simple estimator for $\theta$ is the least squares estimator

$\widehat{\theta}_{T}=\left(\textstyle\sum_{t=1}^{T}\widehat{x}_{t}\widehat{x}_{t-1}^{\mathsf{T}}\right)\left(\textstyle\sum_{t=1}^{T}\widehat{x}_{t-1}\widehat{x}_{t-1}^{\mathsf{T}}\right)^{-1},$
(2)

which may take any value in $\mathbb{R}^{n\times n}$. It is therefore possible
that $\widehat{\theta}_{T}$ is unstable even though $\theta$ is stable, in
which case the estimator is of limited practical value. Alternative estimators
with attractive statistical properties that are guaranteed to be stable have
been proposed in [1, 2, 3, 4]. However, stability is not the only property of
$\theta$ that impacts the qualitative behavior of a linear system; see Figure
1. As optimal control laws are known to inherit important structural
properties from the system matrix [5, Theorem III.1], one should aim to find
estimators that are structurally equivalent to $\theta$. That is, if the least
squares estimator $\widehat{\theta}_{T}$ is structurally different from
$\theta$ itself, in a sense to be made precise later, then implementing
optimal linear feedback designed for $\widehat{\theta}_{T}$ on $\theta$
results in a closed-loop system that is structurally different form the
predicted closed-loop system, e.g., you predict a damper but get a spring.
Capturing the correct qualitative behaviour in the scalar case translates to
enforcing stability and to the need of estimating the sign of $\theta$
correctly. Generalizing this qualitative notion of topological equivalence to
higher dimensions will be the main subject of Section II-A below.

###  I-A Related work

Linear system identification—especially by means of least squares
techniques—has a rich history [6, 7]. In this paper we are, however, not only
interested in finding estimators that fall into the vicinity of the unknown
true model $\theta$. In addition, the estimators should display a
qualitatively similar behavior as $\theta$. This requirement relates to some
extent to the work on qualitative identification pioneered by [8]. More
recently, the focus in linear system identification shifted towards ensuring
the efficient use of data. General informativity of data is discussed in [9],
which justifies the identification pipeline for a class of control problems.
Moreover, sharp statistical characterizations of the effectiveness of the
least squares estimator (2) are presented in [10, 11]. These statistical
results usually quantify the likelihood that $\theta$ lies in some ball around
$\widehat{\theta}_{T}$. However, the models residing within this ball may be
qualitatively different. Leveraging recent results from the theory of large
and moderate deviations [12, 13, 14], we will be able to characterize the
likelihood that the estimated system is qualitatively equivalent to the
unknown true system.

Regarding topological equivalence in the context of linear control systems,
[15] stated in 1980 that “Because of the obvious … practical importance of
these concepts, … there is no doubt that they will become standard vocabulary
among practitioners.” Although [16] later provided many additional insights,
there has been little recent follow-up work on topological properties of
control systems.

###  I-B Contributions

A high-level aim of this work is to showcase how topological insights can
benefit the control community. More specifically, we establish topological
properties of the reverse $I$-projection $\mathcal{P}(\cdot)$ introduced in
[4], which projects any matrix in $\mathbb{R}^{n\times n}$ onto the non-convex
set of stable matrices with respect to an information divergence and can be
evaluated highly efficiently. By exploiting tools from moderate deviations
theory, we characterize here the probability that the reverse $I$-projection
$\mathcal{P}(\widehat{\theta}_{T})$ of the least-squares estimator
$\widehat{\theta}_{T}$ is topologically different from $\theta$. Formally, we
show that

$\mathbb{P}(\mathcal{P}(\widehat{\theta}_{T})\not\overset{t}{\sim}\theta)\lesssim
e^{-\mathcal{O}(\sigma_{\mathrm{min}}(\theta)^{2}\sqrt{T})},$

where ‘$\overset{t}{\sim}$’ denotes topological equivalence. Thus, the
probability that $\mathcal{P}(\widehat{\theta}_{T})$ misrepresents the
topological properties of $\theta$ decays exponentially with $T$ at a rate
$\propto\sigma_{\mathrm{min}}(\theta)^{2}$.

Figure 1: Given any $a\in(0,1)$, the two uncontrolled systems at the top
display structurally different behavior, i.e., they respectively represent a
damper versus a spring. Applying an optimal LQR controller with feedback gain
$k^{\star}$ preserves the structures of the corresponding uncontrolled systems
[5].

Notation. The spectral radius of a matrix $A\in\mathbb{R}^{n\times n}$ is
denoted by $\rho(A)$ and the set of asymptotically stable matrices is denoted
by $\Theta=\\{\theta\in\mathbb{R}^{n\times n}:\rho(\theta)<1\\}$. For a set
$\mathcal{D}\subset\mathbb{R}^{n}$, we use $\mathsf{cl}\,\mathcal{D}$ and
$\mathsf{int}\,\mathcal{D}$ to denote the closure and the interior of
$\mathcal{D}$, respectively. For a real sequence
$\\{a_{T}\\}_{T\in\mathbb{N}}$ the relation $1\ll a_{T}\ll T$ means that
$a_{T}/T\to 0$ and $a_{T}\to\infty$ as $T\to\infty$. We denote the real
$n$-dimensional general linear group by
$\mathsf{GL}(n,\mathbb{R})=\\{A\in\mathbb{R}^{n\times n}:\mathrm{det}(A)\neq
0\\}$. The sets $\mathsf{GL}^{+}(n,\mathbb{R})$ and
$\mathsf{GL}^{-}(n,\mathbb{R})$ contain all matrices in
$\mathsf{GL}(n,\mathbb{R})$ with a strictly positive or negative determinant.

##  II Preliminaries

We first discuss the notion of topological equivalence of linear systems in §
II-A and subsequently review the concept of a moderate deviation principle in
§ II-B.

###  II-A Topological equivalence of linear dynamical systems

This section is mainly inspired by the work of Kuiper and Robbin [17, 18]. In
the following we denote by $f(x)=\theta x$ the state-dependent part of the
dynamics (1), and we refer to the function $f$ as the time-one map. While all
time-one maps considered in this paper are linear, our results naturally
extend to nonlinear systems.

######  Definition II.1 (Topological equivalence).

Two linear time-one maps $f:\mathcal{X}\to\mathcal{X}$ and
$g:\mathcal{Y}\to\mathcal{Y}$ over topological vector spaces $\mathcal{X}$ and
$\mathcal{Y}$ are called topologically equivalent (or conjugate), denoted as
$f\overset{t}{\sim}g$, if there exists a homeomorphism
$\varphi:\mathcal{X}\to\mathcal{Y}$ such that $g\circ\varphi=\varphi\circ f$.

Recall that a homeomorphism is a continuous bijection with a continjous
inverse. We say that two dynamical systems are topologically equivalent if
their time-one maps are topologically equivalent in the sense of Definition
II.1. If the two dynamical systems are noise-free (i.e., the time-one maps
alone determine the dynamics), then they are topologically equivalent if there
is a homeomorphism mapping the trajectories of one system onto those of the
other [19], [20, Chapter 2]. This means that the state trajectories of
topologically equivalent dynamical systems are qualitatively identical. To see
this, consider two time-one maps $f:\mathcal{X}\to\mathcal{X}$ and
$g:\mathcal{Y}\to\mathcal{Y}$ initialized at $x_{0}$ and $y_{0}$,
respectively, and assume that $f$ and $g$ are topologically equivalent. In
this case

${x_{0}}$${f^{n}(x_{0})}$${y_{0}}$${g^{n}(y_{0})}$$\scriptstyle{f^{n}}$$\scriptstyle{\varphi}$$\scriptstyle{\varphi}$$\scriptstyle{g^{n}}$

commutes for any $n\in\mathbb{N}$ because $f=\varphi\circ g\circ\varphi^{-1}$
[20]. Note that this argument holds indeed for any pair $(x_{0},y_{0})$ of
initial states, that is, $\varphi$ essentially constitutes a homeomorphic
coordinate change. A change of a dynamical system that destroys topological
equivalence is called a bifurcation. If two linear time-one maps of the form
$f(x)=Fx$ and $g(y)=Gy$ parametrized by the matrices $F$ and $G$ are
topologically equivalent, then, by slight abuse of notation, we write
$F\overset{t}{\sim}G$. We emphasize that topological equivalence generalizes
the more common _linear equivalence_ , where $\varphi(x)=Zx$ is a linear
isomorphism parametrized by a matrix $Z$. In this case the time-one maps
$f(x)=Fx$ and $g(y)=Gy$ with $G=ZAZ^{-1}$ are linearly equivalent. To gain
some intuition for Definition II.1, consider a scalar system of the form
$f(x)=ax$ for some $a\in\mathbb{R}$. Figure 2 shows that this system admits
seven equivalence classes with respect to the topological equivalence relation
at hand, which correspond to seven intervals in $\mathbb{R}$. This example
also shows that one should not expect both $\varphi$ and $\varphi^{-1}$ to be
differentiable. Indeed, assuming $\varphi$ to be a diffeomorphism implies that
$\varphi$ must be linear [19, Proposition 6.1]. Note also that any
homeomorphism on $\mathbb{R}$ is necessarily monotone, which confirms our
earlier insight that the damper can never be mapped to the spring in Figure 1.

Figure 2: Visualization of the seven topological equivalence classes of a
scalar linear system [18, Proposition 1.5]. Note that if the time-one maps
$f(x)=ax$ and $g(y)=by$ belong to the same non-degenerate class (i.e., $a$ and
$b$ both belong to interval $(1)$, $(3)$, $(5)$ or $(7)$), then
$g=\varphi\circ f\circ\varphi^{-1}$ for $\varphi(x)=x|x|^{c-1}$ with
$c=\log(|b|)/\log(|a|)$. For example, $f(x)=2x$ and $g(y)=8y$ belong to class
$(7)$ and are related by the smooth map $\varphi(x)=x^{3}$ with inverse
$\varphi^{-1}=x^{1/3}$ [18].

Asymptotic stability or topological invariants such as orientation provide
conditions for the topological equivalence of two time-one maps [18].

######  Definition II.2 (Orientation).

We call a linear invertible time-one map $f(x)=Fx$ orientation-preserving if
$\mathrm{det}(F)>0$, and we denote by $\mathrm{or}(f)$ the sign of
$\mathrm{det}(F)$.

Note that $\mathsf{GL}^{+}(n,\mathbb{R})$ and $\mathsf{GL}^{-}(n,\mathbb{R})$
represent the components of $\mathsf{GL}(n,\mathbb{R})$ that contain all
invertible matrices $F$ whose associated time-one maps $f(x)=Fx$ satisfy
$\mathrm{or}(f)=1$ and $\mathrm{or}(f)=-1$, respectively. Note also that the
sign of the (signed) volume of the unit hypercube is invariant under an
orientation-preserving map. For additional information on the concept of
orientation see for example [21, Chapter 6]. The orientation of a map is also
known to be a _topological invariant_ , that is, if two invertible maps $f$
and $g$ are topologically equivalent, then $\mathrm{or}(f)=\mathrm{or}(g)$
[22, Chapters 6 & 10]. For brevity, we do not generalize the orientation
operator to linear maps that fail to be invertible. The following lemma
provides a convenient tool for checking topological equivalence.

######  Lemma II.3 (Topological equivalence of asymptotically stable systems
[19, Theorem 9.2] ).

Assume that $f(x)=Fx$ and $g(y)=Gy$ are asymptotically stable linear
isomorphisms on $\mathbb{R}^{n}$. Then, we have $f\overset{t}{\sim}g$ if and
only if there exists a continuous path $X(t)$ in $\mathsf{GL}(n,\mathbb{R})$
parametrized in $t\in[0,1]$ with $X(0)=F$ and $X(1)=G$.

As $\mathsf{GL}(n,\mathbb{R})$ can be decomposed into
$\mathsf{GL}^{+}(n,\mathbb{R})$ and $\mathsf{GL}^{-}(n,\mathbb{R})$, both of
which are connected, Lemma II.3 implies that these $f$ and $g$ are
topologically equivalent if and only if the signs of the determinants of $F$
and $G$ match. Be slight abuse of notation, we henceforth define the
orientation $\mathrm{or}(F)$ of an invertible matrix $F$ as the sign of
$\det(F)$.

###  II-B Moderate deviations theory

Throughout the paper we assume that all random objects are defined on a
measurable space $(\Omega,\mathcal{F})$ equipped with a probability measure
$\mathbb{P}_{\theta}$ parametrized by the (asymptotically stable) system
matrix $\theta\in\Theta$ from (1). We denote the expectation operator with
respect to $\mathbb{P}_{\theta}$ by $\mathbb{E}_{\theta}[\cdot]$. From now on
we impose the following assumption borrowed from [4].

######  Assumption II.4 (Linear system).

The following hold.

  1. (i)

The system (1) is asymptotically stable, i.e., $\theta\in\Theta$.

  2. (ii)

For each $\theta\in\Theta$ the disturbances $\\{w_{t}\\}_{t\in\mathbb{N}}$ are
independent and identically distributed (i.i.d.) and independent of $x_{0}$
under $\mathbb{P}_{\theta}$. The disturbances are unbiased
($\mathbb{E}_{\theta}[w_{t}]=0$) and non-degenerate
($S_{w}=\mathbb{E}_{\theta}[w_{t}w_{t}^{\mathsf{T}}]\succ 0$), and their
probability density is everywhere positive.

Assumption II.4 implies that the linear system (1) admits an invariant
distribution $\nu_{\theta}$ [23, § 10.5.4] with zero mean and covariance
matrix $S_{\theta}$, which is given by the unique positive definite solution
of the discrete Lyapunov equation

$S_{\theta}=\theta S_{\theta}\theta^{\mathsf{T}}+S_{w},$ (3)

see, e.g., [24, § 6.10 E]. Next, we describe a method to characterize the
probability of the least squares estimator $\widehat{\theta}_{T}$ deviating
from $\theta$ by a prescribed threshold. To this end, we denote by
$\Theta^{\prime}=\mathbb{R}^{n\times n}$ the space of all estimator
realizations that are possible in view of Assumption II.4, and we use the
discrepancy function $I:\Theta^{\prime}\times\Theta\rightarrow[0,\infty]$ with

$I(\theta^{\prime},\theta)=\frac{1}{2}\mathrm{tr}\left(S_{w}^{-1}(\theta^{\prime}-\theta)S_{\theta}(\theta^{\prime}-\theta)^{\mathsf{T}}\right)$
(4)

to quantify the difference between an estimator realization
$\theta^{\prime}\in\Theta^{\prime}$ and the system matrix $\theta\in\Theta$;
see [4]. Here, the invariant state covariance matrix $S_{\theta}$ is defined
as in (3). Note that $S_{\theta}$ and thus also $I(\theta^{\prime},\theta)$
diverge as $\theta$ approaches the boundary of $\Theta$ and becomes unstable.
Note also that since $S_{w}\succ 0$ and hence $S_{\theta}\succ 0$,
$I(\theta^{\prime},\theta)$ vanishes if and only if $\theta^{\prime}=\theta$.
In this sense $I$ behaves like a distance. Note, however, that
$I(\theta^{\prime},\theta)$ is not symmetric in $\theta$ and
$\theta^{\prime}$.

Next we recall the notions of a rate function and a moderate deviation
principle (MDP) such that we can review the key results of [4]. For a
comprehensive introduction to moderate deviations theory we refer to [13, 14].

######  Definition II.5 (Rate function).

We call $I:\Theta^{\prime}\times\Theta\rightarrow[0,\infty]$ a rate function
if it is lower semi-continuous in $\theta^{\prime}$.

######  Definition II.6 (Moderate deviation principle).

A sequence of estimators $\\{\widehat{\vartheta}_{T}\\}_{T\in\mathbb{N}}$ is
said to satisfy a moderate deviation principle with rate function $I$ if for
every sequence $\\{a_{T}\\}_{T\in\mathbb{N}}$ of real numbers with $1\ll
a_{T}\ll T$, for every Borel set $\mathcal{D}\subset\Theta^{\prime}$ and for
every $\theta\in\Theta$ we have

$\displaystyle-\inf_{\theta^{\prime}\in\mathsf{int}\,{\mathcal{D}}}\,I(\theta^{\prime},\theta)\leq$
$\displaystyle\liminf_{T\to\infty}~{}\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\left(\widehat{\vartheta}_{T}\in\mathcal{D}\right)$
(5a) $\displaystyle\leq$
$\displaystyle\limsup_{T\to\infty}~{}\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\left(\widehat{\vartheta}_{T}\in\mathcal{D}\right)$
(5b) $\displaystyle\leq$
$\displaystyle-\inf_{\theta^{\prime}\in\mathsf{cl}\,{\mathcal{D}}}\,I(\theta^{\prime},\theta).$
(5c)

If $I(\theta^{\prime},\theta)$ is continuous (instead of merely lower semi-
continuous) in $\theta^{\prime}$ and if $\mathsf{int}\,\mathcal{D}$ is dense
in $\mathcal{D}$, then the infima in (5a) and (5c) match, and all inequalities
in (5) reduce to equalities. Then, (5) implies that
$\mathbb{P}_{\theta}(\widehat{\vartheta}_{T}\in\mathcal{D})=e^{-ra_{T}+o(a_{T})}$,
where $r=\inf_{\theta^{\prime}\in\mathcal{D}}I(\theta^{\prime},\theta)$
quantifies the $I$-distance of the system matrix $\theta$ from the set
$\mathcal{D}$, see Figure 3. Thus, the distance $r$ coincides with the decay
rate of the probability that $\widehat{\vartheta}_{T}$ materializes in
$\mathcal{D}$, while $\\{a_{T}\\}_{T\in\mathbb{N}}$ can be viewed as the speed
of convergence. The condition $1\ll a_{T}\ll T$ is satisfied, for example, if
$a_{T}=\sqrt{T}$, $T\in\mathbb{N}$.

Figure 3: The $I$-distance between $\theta$ and some set
$\mathcal{D}\not\ni\theta$ quantifies the rate at which the likelihood of an
estimator $\widehat{\theta}_{T}$ falling into $\mathcal{D}$ decays, i.e.,
$\mathbb{P}_{\theta}(\widehat{\theta}_{T}\in\mathcal{D})\approx e^{-ra_{T}}$.

It is unknown whether the least squares estimators $\widehat{\theta}_{T}$
defined in (2) satisfy an MDP for $n>1$. However, the transformed least
squares estimators defined as

$\widehat{\vartheta}_{T}(\widehat{\theta}_{T},\theta)=\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)+\theta$
(6)

are known to satisfy an MDP for any $n\in\mathbb{N}$, where the underlying
rate function given by the discrepancy function (4). Whenever there is no risk
of confusion, we drop the explicit dependence of $\widehat{\vartheta}_{T}$ on
$\widehat{\theta}_{T}$ and $\theta$. In order to formally introduce the
advertised MDP, we impose another standard regularity condition borrowed from
[4], which is again assumed to hold throughout the remainder of the paper.

######  Assumption II.7 (Light-tailed noise and stationarity).

The following hold for every $\theta\in\Theta$.

  1. (i)

The disturbances $\\{w_{t}\\}_{t\in\mathbb{N}}$ are light-tailed, i.e., there
exists $\alpha>0$ with $\mathbb{E}_{\theta}[e^{\alpha\|w_{t}\|^{2}}]<\infty$
for all $t\in\mathbb{N}$.

  2. (ii)

The initial distribution $\nu$ coincides with the invariant distribution
$\nu_{\theta}$ of the linear system (1).

We can now formally state the MDP for the estimators (6).

######  Proposition II.8 (Moderate deviation principle [4, Proposition 3.4]).

If $\\{a_{T}\\}_{T\in\mathbb{N}}$ is a real sequence with $1\ll a_{T}\ll T$,
then the transformed least squares estimators
$\\{\widehat{\vartheta}_{T}\\}_{T\in\mathbb{N}}$ defined in (6) satisfy an MDP
with rate function (4).

Proposition II.8 is attractive due to the generality of MDPs, which provide
tight bounds on the (asymptotic) probability of any Borel set of estimator
realizations. A simple direct application of Proposition II.8 is described
below.

######  Example II.9 (Scalar system identification and noise invariance).

Consider a scalar system with $S_{w}=\sigma_{w}^{2}>0$. As shown in [4,
Example 3.5], Proposition II.8 implies that

$\displaystyle\mathbb{P}_{\theta}(|\widehat{\theta}_{T}-\theta|$
$\displaystyle\\!>\varepsilon\sqrt{a_{T}/T})$
$\displaystyle=\textstyle\exp\left(-\frac{1}{2}\varepsilon^{2}\,a_{T}/(1-\theta^{2})+o(a_{T})\right)$

for any $\varepsilon>0$ and $T\in\mathbb{N}$. Thus, the decay rate on right
hand side of the above expression is independent of the noise intensity
$\sigma_{w}^{2}$. If $n>1$, define
$Z(\theta^{\prime},\theta)=(\theta^{\prime}-\theta)\otimes(\theta^{\prime}-\theta)(I_{n^{2}}-\theta\otimes\theta)^{-1}$,
where $\otimes$ denotes the Kronecker product. By [25, p. 265], the rate
function (4) can then be recast as

$\displaystyle I(\theta^{\prime},\theta)$
$\displaystyle=\frac{1}{2}\mathrm{vec}(S_{w}^{-1})^{\mathsf{T}}Z(\theta^{\prime},\theta)\mathrm{vec}(S_{w}),$

where $\mathrm{vec}(S_{w})\in\mathbb{R}^{n^{2}}$ represents the vector
obtained by stacking the columns of $S_{w}$ on top of each other. This
representation reveals that the rate function is indeed invariant under
scaling of the noise covariance matrix, that is, it is independent of the
overall noise level for all $n\in\mathbb{N}$. ∎

We finally highlight that the rate function (4) can be used to construct a
reverse $I$-projection defined through

$\mathcal{P}(\theta^{\prime})\in\arg\min_{\theta\in\Theta}I(\theta^{\prime},{\theta}),$
(7)

which maps any point $\theta^{\prime}\in\Theta^{\prime}$ to the nearest point
in the non-convex set $\Theta$ of asymptotically stable matrices with respect
to the $I$-distance; see [4]. The reverse $I$-projection
$\mathcal{P}(\widehat{\theta}_{T})$ of the least squares estimator
$\widehat{\theta}_{T}$ provides an estimator for $\theta$ that is stable by
construction, has desirable statistical properties and can be efficiently
computed. Indeed, one can show that

$\mathcal{P}(\theta^{\prime})=\theta^{\prime}+\mathsf{dlqr}(\theta^{\prime},I_{n},I_{n},(2\delta
S_{w})^{-1})+\mathcal{O}(\delta^{p})$ (8)

for some $p\geq 1$, where $\mathsf{dlqr}(\cdot)$ denotes the standard
discrete-time LQR
routine,111https://juliacontrol.github.io/ControlSystems.jl/latest/examples/example/#LQR-
design see [4, §3.3]. In addition, the reverse $I$-projection preserves
orientation, i.e.,
$\mathrm{or}(\mathcal{P}(\theta^{\prime}))=\mathrm{or}(\theta^{\prime})$ for
any $\theta^{\prime}\in\mathsf{GL}(n,\mathbb{R})$, see, e.g., [4, Corollary
3.12] and [5]. In fact, the numerical approximation of
$\mathcal{P}(\theta^{\prime})$ on the right hand side of (8) without the error
term $\mathcal{O}(\delta^{p})$ is also asymptotically stable and preserves
orientation for any $\delta>0$. Due to its desirable statistical and
computational properties, our proposed approach to estimate the topological
class of $\theta\in\Theta$ will critically rely on the reverse $I$-projection
$\mathcal{P}(\widehat{\theta}_{T})$.

##  III Main results

To aid the presentation we assume from now on without much loss of generality
that $\theta$ is invertible. We are now ready to demonstrate that the MDP of
Proposition II.8 allows us via the reverse $I$-projection
$\mathcal{P}(\widehat{\theta}_{T})$ to derive sharp bounds on the decay rate
of the probability of the event
$\mathcal{P}(\widehat{\theta}_{T})\not\overset{t}{\sim}\theta$. Recall from
Lemma II.3 that the two stable and ($\mathbb{P}_{\theta}$-almost surely)
invertible matrices $\mathcal{P}(\widehat{\theta}_{T})$ and $\theta$ are
topologically equivalent if and only if they have the same orientation. Recall
also that
$\mathrm{or}(\mathcal{P}(\widehat{\theta}_{T}))=\mathrm{or}(\widehat{\theta}_{T})$
because the reverse $I$-projection preserves orientation. Checking whether
$\mathcal{P}(\widehat{\theta}_{T})$ is topologically equivalent to $\theta$ is
thus tantamount to checking whether the determinants of $\widehat{\theta}_{T}$
and $\theta$ have the same signs.

######  Theorem III.1 (Probability of misclassification).

Assume that $\theta\in\Theta\cap\mathsf{GL}(n,\mathbb{R})$,
$\\{\widehat{\theta}_{T}\\}_{T\in\mathbb{N}}$ are the least squares estimators
(2) and $\\{a_{T}\\}_{T\in\mathbb{N}}$ is a sequence with $1\ll a_{T}\ll T$.
If $S_{\theta^{\circ}}=S_{w}^{-1/2}S_{\theta}S_{w}^{-1/2}$ and
$r=\frac{1}{2}\lambda_{\mathrm{min}}(S_{\theta^{\circ}}-I_{n})$, then

$\limsup_{T\to\infty}\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\big{(}\mathrm{or}(\widehat{\theta}_{T})\not=\mathrm{or}(\theta)\big{)}\leq-r,$
(9a)
$\limsup_{T\to\infty}\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\big{(}\mathcal{P}(\widehat{\theta}_{T})\not\overset{t}{\sim}\theta\big{)}\leq-r.$
(9b)

###### Proof.

As for (9a), note that $\widehat{\theta}_{T}\in\mathsf{GL}(n,\mathbb{R})$
$\mathbb{P}_{\theta}$-almost surely for all sufficiently large $T$, and
therefore we have

$\displaystyle\mathbb{P}_{\theta}\big{(}\mathrm{or}(\widehat{\theta}_{T})\not=\mathrm{or}(\theta)\big{)}$
$\displaystyle=\mathbb{P}_{\theta}\big{(}\mathrm{or}(\widehat{\theta}_{T})=-\mathrm{or}(\theta)\big{)}$
$\displaystyle=\mathbb{P}_{\theta}\big{(}\exists
G\in\mathsf{GL}^{-}(n,\mathbb{R}):\widehat{\theta}_{T}=G\theta\big{)}$
$\displaystyle=\mathbb{P}_{\theta}\big{(}\widehat{\vartheta}_{T}\in\mathcal{D}_{T}(\theta)\big{)}$

where the second equality holds because any invertible matrices
$\widehat{\theta}_{T}$ and $\theta$ whose determinants have opposite signs
satisfy $\widehat{\theta}_{T}=G\theta$ for some
$G\in\mathsf{GL}^{-}(n,\mathbb{R})$. The third equality follows from the
definition of the transformed least squares estimators
$\\{\widehat{\vartheta}_{T}\\}_{T\in\mathbb{N}}$ in (6) and the construction
of the set

$\displaystyle\mathcal{D}_{T}(\theta)=\\{$
$\displaystyle\sqrt{(T/a_{T})}(G-I_{n})\theta+\theta:G\in\mathsf{GL}^{-}(n,\mathbb{R})\\}.$

As this set is time-dependent, we cannot directly use it. That is, it is not
admissible in the sense of Definition II.6. To sidestep this complication, we
consider instead the larger set

$\mathcal{D}(\theta)=\bigcup_{T\in\mathbb{N}}\mathcal{D}_{T}(\theta).$

The above reasoning then implies that

$\displaystyle\limsup_{T\to\infty}\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\big{(}\mathrm{or}(\widehat{\theta}_{T})\not=\mathrm{or}(\theta)\big{)}$
$\displaystyle\leq$
$\displaystyle\limsup_{T\to\infty}\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\big{(}\widehat{\vartheta}_{T}\in\mathcal{D}(\theta)\big{)}\leq-\inf_{\theta^{\prime}\in\mathsf{cl}\,\mathcal{D}(\theta)}I(\theta^{\prime},\theta),$

where the first inquality holds because
$\mathcal{D}_{T}(\theta)\subset\mathcal{D}(\theta)$ for all $T\in\mathbb{N}$,
while the second inequality follows form the MDP established in Proposition
II.8. In the remainder of the proof we derive an analytical lower bound on the
minimization problem on the right hand side of the above expression. To this
end, assume first that $S_{w}=I_{n}$. Evaluating the rate function (4) at an
arbitrary $\theta^{\prime}\in\mathcal{D}(\theta)$ then yields

$\displaystyle
I(\theta^{\prime},\theta)=\frac{T}{2a_{T}}\mathsf{tr}\left((G-I_{n})\theta
S_{\theta}\theta^{\mathsf{T}}(G-I_{n})^{\mathsf{T}}\right)$ (10)

for some $G\in\mathsf{GL}^{-}(n,\mathbb{R})$ and some $T\in\mathbb{N}$, and
the Lyapunov equation (3) implies that $\theta
S_{\theta}\theta^{\mathsf{T}}=S_{\theta}-I_{n}$. In addition, our assumptions
about the sequence $\\{a_{T}\\}_{T\in\mathbb{N}}$ imply that $T/a_{T}\geq 1$.
We may thus conclude that

$\displaystyle\min_{\theta^{\prime}\in\mathsf{cl}\,\mathcal{D}(\theta)}$
$\displaystyle I(\theta^{\prime},\theta)$ $\displaystyle\geq\inf_{\det(G)\leq
0}\frac{1}{2}\mathsf{tr}\left((G-I_{n})(S_{\theta}-I_{n})(G-I_{n})^{\mathsf{T}}\right)$
$\displaystyle\geq\frac{1}{2}(\lambda_{\mathrm{min}}(S_{\theta})-1)\inf_{\det(G)\leq
0}\|G-I_{n}\|_{\mathsf{F}}^{2}$
$\displaystyle=\frac{1}{2}(\lambda_{\mathrm{min}}(S_{\theta})-1),$

where the first inequality holds because $\det(G)\leq 0$ for every
$G\in\mathsf{cl}\,\mathsf{GL}^{-}(n,\mathbb{R})$, the second inequality uses
the bound $\mathsf{tr}(AB)\geq\sigma_{\mathrm{min}}(A)\mathsf{tr}(B)$ for any
$A,B\succeq 0$, and the third inequality follows from the Eckart-Young Theorem
[26, Theorem 2.4.8]. One can actually show that the second inequality is
tight, but this is not needed to prove the theorem. This establishes (9a) for
$S_{w}=I_{n}$. If $S_{w}\succ 0$ is arbitrary, one may first apply a change of
coordinates $x^{\circ}=S_{w}^{-1/2}x$, under which the noise covariance matrix
simplifies to $I_{n}$, while the system matrix and the invariant state
covariance matrix become $\theta^{\circ}=S_{w}^{-1/2}\theta S_{w}^{1/2}$ and
$S_{\theta^{\circ}}=S_{w}^{-1/2}S_{\theta}S_{w}^{-1/2}$, respectively. As
$S_{w}\succ 0$, we further have $\theta^{\circ}\overset{t}{\sim}\theta$.
Applying the results of the first part of the proof to the transformed system
finally yields (9a). As for (9b), recall from Sections II-A and II-B that
$\mathcal{P}(\theta^{\prime})$ is asymptotically stable and that
$\mathrm{or}(\theta^{\prime})=\mathrm{or}(\mathcal{P}(\theta^{\prime}))$ for
any $\theta^{\prime}\in\mathsf{GL}(n,\mathbb{R})$. Therefore, we can focus on
orientation solely, i.e., if
$\mathrm{or}(\theta^{\prime})=\mathrm{or}(\theta)$ then
$\mathcal{P}(\theta^{\prime})\overset{t}{\sim}\theta$. Then again, as
$\mathbb{P}_{\theta}\big{(}\widehat{\theta}_{T}\not\in\mathsf{GL}(n,\mathbb{R})\big{)}=0$
for sufficiently large $T$, the claim follows from (9a). ∎

The rate $r=\frac{1}{2}\lambda_{\mathrm{min}}(S_{\theta^{\circ}}-I_{n})$
established in Theorem III.1 is non-trivial (strictly positive) for any
$\theta\in\Theta\cap\mathsf{GL}(n,\mathbb{R})$ as

$S_{\theta^{\circ}}=S_{w}^{-1/2}S_{\theta}S_{w}^{-1/2}=S_{w}^{-1/2}\sum^{\infty}_{k=1}\theta^{k}S_{w}(\theta^{k})^{\mathsf{T}}S_{w}^{-1/2}+I_{n},$

where the second equality holds because $S_{\theta}$ satisfies (3).
Quantitative notions of stability and controllability, which strengthen the
standard qualitative notions of stability and controllability, respectively,
offer further insights into (9b).

######  Definition III.2 (Strong controllability [27]).

A pair $(A,B)$ is called $(\ell,\nu)$-strongly controllable for
$\ell\in\mathbb{N}$ and $\nu>0$ if the matrix
$C_{\ell}=\begin{pmatrix}B\,AB\,\cdots\,A^{\ell-1}B\end{pmatrix}$ satisfies
$\sigma_{\mathrm{min}}(C_{\ell})\geq\nu$.

If $A$ and $B$ parametrize the deterministic system $x_{t+1}=Ax_{t}+Bu_{t}$
and $(A,B)$ is $(\ell,\nu)$-strongly controllable, then, a small $\nu$
indicates that for two points $x_{a},x_{b}$ being close, the inputs to reach
them from $x_{0}$ might be wildly different.

######  Definition III.3 (Strong stability [27]).

A square matrix $A$ is $(\kappa,\gamma)$-strongly stable for $\kappa>0$ and
$\gamma\in(0,1]$ if $A=HLH^{-1}$ for some matrices $L$ and $H$ with
$\|L\|_{2}\leq 1-\gamma$ and $\|H\|_{2}\|H^{-1}\|_{2}\leq\kappa$.

One readily verifies that any $(\kappa,\gamma)$-strongly stable matrix is
asymptotically stable. Conversely, by Lyapunov stability theory, for any
asymptotically stable matrix $A$ there exists $P\succ 0$ such that
$P-A^{\mathsf{T}}PA\succ 0$. Thus, the transformed matrix
$A^{\prime}=P^{1/2}AP^{-1/2}$ satisfies
$I_{n}-(A^{\prime})^{\mathsf{T}}A^{\prime}\succ 0$, which in turn implies that
$\|A^{\prime}\|_{2}<1$. Hence, $A=HLH^{-1}$ with $H=P^{-1/2}$ and
$L=A^{\prime}$, and $A$ is $(\kappa,\gamma)$-strongly stable with
$\gamma=1-\|A^{\prime}\|_{2}$ and $\kappa$ being proportional to the condition
number of $P^{1/2}$. The next lemma relates the rate $r$ of Theorem III.1 to
strong controllability and stability properties of the pair
$(\theta^{\circ},I_{n})$.

######  Lemma III.4 (Bounding the rate of Theorem III.1).

If the pair $(\theta^{\circ},I_{n})$ is $(\ell,\nu)$-strongly controllable and
$\theta^{\circ}$ is $(\kappa,\gamma)$-strongly stable, then

$\nu^{2}\leq\lambda_{\mathrm{min}}(S_{\theta^{\circ}})\leq\frac{\kappa^{2}}{2\gamma-\gamma^{2}}.$
(11)

###### Proof.

The matrix $C_{\ell}$ introduced in Definition III.2 with $A=\theta^{\circ}$
and $B=I_{n}$ satisfies
$C_{\ell}C_{\ell}^{\mathsf{T}}=\sum^{\ell-1}_{i=0}(\theta^{\circ})^{i}((\theta^{\circ})^{i})^{\mathsf{T}}$.
Thus, the controllability Gramian
$W=\lim_{\ell\to\infty}C_{\ell}C_{\ell}^{\mathsf{T}}$ coincides with
$S_{\theta^{\circ}}$ as defined in Theorem III.1. As $(\theta^{\circ},I_{n})$
is $(\ell,\nu)$-strongly controllable, we may thus conclude that
$\lambda_{\mathrm{min}}(S_{\theta^{\circ}})\geq\lambda_{\mathrm{min}}(C_{\ell}C_{\ell}^{\mathsf{T}})\geq\nu^{2}$.
An upper bound on $\lambda_{\mathrm{min}}(S_{\theta^{\circ}})$ can be obtained
from the $(\kappa,\gamma)$-strong stability of $\theta^{\circ}$, which implies
via [27, Lemma 3.3] that
$\mathsf{tr}(S_{\theta^{\circ}})\leq(\kappa^{2}/\gamma)\mathsf{tr}(I_{n})$ and
hence $\lambda_{\mathrm{min}}(S_{\theta^{\circ}})\leq(\kappa^{2}/\gamma)$.
Looking at the proof of [27, Lemma 3.3], this can be sharpened to (11). ∎

The bounds (11) confirm the folklore wisdom in system identification that a
slow ($\gamma$ small), yet well-excited ($\nu$ large) system is desirable.
Next, we highlight a few other properties of the problem parameters
$(\theta,S_{w})$ one might be able to manipulate in order to increase the rate
in Theorem III.1.

###  III-A Tuning $\sigma_{\mathrm{min}}(\theta^{\circ})$ and the noise
covariance matrix

Using the series representation of $S_{\theta^{\circ}}$ one can show that

$S_{\theta^{\circ}}\succeq
I_{n}+\frac{1}{1-\sigma_{\mathrm{min}}(\theta^{\circ})^{2}}\theta^{\circ}(\theta^{\circ})^{\mathsf{T}},$

see [28], where
$\sigma_{\mathrm{min}}(\theta^{\circ})\leq\rho(\theta^{\circ})<1$. Since also
$S_{\theta^{\circ}}-I_{n}\succeq\theta^{\circ}(\theta^{\circ})^{\top}$ it
follows that

$\frac{1}{2}\lambda_{\mathrm{min}}(S_{\theta^{\circ}}-I_{n})\geq\frac{1}{2}\sigma_{\mathrm{min}}(\theta^{\circ})^{2}.$

Therefore, an increase in $\sigma_{\mathrm{min}}(\theta^{\circ})$ thus
improves the rate $r$ from Theorem III.1. Next, assume that $\theta$ is
diagonalizable, that is, $\theta=V\Lambda V^{-1}$ for some diagonal matrix
$\Lambda$ and invertible matrix $V$. As
$\sigma_{\mathrm{min}}(\theta^{\circ})\leq\lambda_{\mathrm{min}}(\theta^{\circ})=\lambda_{\mathrm{min}}(\theta)$,
the preferred noise covariance matrix for which
$\sigma_{\mathrm{min}}(\theta^{\circ})$ matches the bound
$\lambda_{\mathrm{min}}(\theta)$ (which is independent of $S_{w}$) is given by
$S_{w}=\alpha VV^{\mathsf{T}}$ for any $\alpha>0$. As already pointed out in
Example II.9 the overall magnitude of $S_{w}$ is not important, but its
principal axes are.

###  III-B Tuning $\sigma_{\mathrm{min}}(\theta)$

The rate $\mathcal{O}(-\sigma_{\mathrm{min}}(\theta)^{2}\sqrt{T})$ from the
introduction follows as
$\sigma_{\mathrm{min}}(\theta)\leq\sigma_{\mathrm{min}}(\theta^{\circ})\sigma_{\mathrm{max}}(S_{w}^{1/2})/\sigma_{\mathrm{min}}(S_{w}^{1/2})$.
In a system-theoretic context, as
$\sigma_{\mathrm{min}}(\theta)=\inf_{\|x\|_{2}=1}\|\theta x\|_{2}$, an
increase in $\sigma_{\mathrm{min}}(\theta)$ reduces the contraction rate of
(1).

###  III-C Minimizing interconnections

Consider a separable and an interconnected system with

$\theta_{1}=\begin{pmatrix}Y&0\\\
0&Y\end{pmatrix}\quad\text{and}\quad\theta_{2}=\begin{pmatrix}Y&I_{n}\\\
0&Y\end{pmatrix}$

for some $Y\in\Theta$, respectively. As

$\theta_{1}^{\mathsf{T}}\theta_{1}=\begin{pmatrix}Y^{\mathsf{T}}Y&0\\\
0&Y^{\mathsf{T}}Y\end{pmatrix}~{}\text{and}~{}\theta_{2}^{\mathsf{T}}\theta_{2}=\begin{pmatrix}Y^{\mathsf{T}}Y&Y^{\mathsf{T}}\\\
Y&Y^{\mathsf{T}}Y\end{pmatrix},$

one can show that
$\sigma_{\mathrm{min}}(\theta_{1})\geq\sigma_{\mathrm{min}}(\theta_{2})$. This
aligns with intuition: if possible, identify subsystems individually.

##  IV Numerical Example

We now compare the theoretical decay rate of topological misclassification
derived in Theorem III.1 against the empirical decay rate for the nominal
least squares estimator $\widehat{\theta}_{T}$ and its reverse $I$-projection
$\mathcal{P}(\widehat{\theta}_{T})$. Concurrently, we exemplify the insights
from Section III-C. To this end, we set

$Y=\begin{pmatrix}-0.1&1\\\ 0.1&0.05\end{pmatrix},$

and simulate system (1) for both $\theta_{1}$ and $\theta_{2}$ defined as in
Section III-C starting from $E=10^{3}$ initial conditions
$x_{0}\overset{i.i.d.}{\sim}\mathcal{N}(0,I_{4})$ with $S_{w}=I_{4}$. Each
initial condition leads to a trajectory under both $\theta_{1}$ and
$\theta_{2}$ from which we construct the corresponding least squares
estimators $\widehat{\theta}^{(i)}_{j,T}$, $i=1,\dots,E$, $T=1,\dots,10^{3}$,
$j\in\\{1,2\\}$. Averaging over the $E$ simulation runs yields the empirical
probability that $\widehat{\theta}_{j,T}$ or its reverse $I$-projection are
topologically equivalent to the true system matrix $\theta_{j}$.

Figure 4: Empirical versus theoretical convergence rates.

Figure (4) compares the bounds on the misclassification probability derived in
Theorem III.1 for $a_{T}=T^{\frac{1}{1+\epsilon}}$ with $\epsilon=10^{-9}$
against the empirical probabilities. Here,
$\mathcal{P}(\widehat{\theta}_{j,T}^{(i)})$ is computed via (8) for
$\delta=10^{-9}$. As expected, the projection accelerates topological
identification, and block-diagonal system matrices are easier to identify. All
experiments are implemented in MATLAB.

## References

  * [1] J.M. Maciejowski  “Guaranteed stability with subspace methods”  In _Systems & Control Letters_ 26.2, 1995, pp. 153–156
  * [2] Byron Boots, Geoffrey J Gordon and Sajid M. Siddiqi  “A Constraint Generation Approach to Learning Stable Linear Dynamical Systems”  In _Advances in Neural Information Processing Systems_ , 2008, pp. 1329–1336
  * [3] Nicolas Gillis, Michael Karow and Punit Sharma  “Approximating the nearest stable discrete-time system”  In _Linear Algebra and its Applications_ 573, 2019, pp. 37–53
  * [4] Wouter Jongeneel, Tobias Sutter and Daniel Kuhn  “Efficient learning of a linear dynamical system with stability guarantees”, 2021  arXiv:2102.03664
  * [5] W. Jongeneel and D. Kuhn  “On Topological Equivalence in Linear Quadratic Optimal Control”  In _submitted_ , 2020
  * [6] Peter Van Overschee and Bart De Moor  “Subspace Identification for Linear Systems: Theory – Implementation – Applications”  Kluwer Academic Publishers, 1996
  * [7] Michel Verhaegen and Vincent Verdult  “Filtering and System Identification”  Cambridge University Press, 2007
  * [8] Benjamin Kuipers  “Qualitative Reasoning: Modeling and Simulation with Incomplete Knowledge”  MIT Press, 1994
  * [9] Henk J Van Waarde, Jaap Eising, Harry L Trentelman and M Kanat Camlibel  “Data informativity: a new perspective on data-driven analysis and control”  In _IEEE Transactions on Automatic Control_ 65.11 IEEE, 2020, pp. 4753–4768
  * [10] Max Simchowitz et al.  “Learning Without Mixing: Towards A Sharp Analysis of Linear System Identification”  In _Conference On Learning Theory_ , 2018, pp. 439–473
  * [11] Tuhin Sarkar and Alexander Rakhlin  “Near optimal finite time identification of arbitrary linear dynamical systems”  In _International Conference on Machine Learning_ , 2019, pp. 5610–5618
  * [12] Paul Dupuis and Richard S. Ellis  “A Weak Convergence Approach to the Theory of Large Deviations”  John Wiley & Sons, 1997
  * [13] F. Hollander  “Large Deviations”  American Mathematical Society, 2008
  * [14] A. Dembo and O. Zeitouni  “Large Deviations Techniques and Applications”  Springer, 2009
  * [15] Jan C Willems  “Topological classification and structural stability of linear systems”  In _Journal of Differential Equations_ 35.3, 1980, pp. 306–318
  * [16] Jan Willem Polderman  “Adaptive Control & Identification: Conflict or Conflux”  PhD thesis, 1987
  * [17] J. W. Robbin  “Topological conjugacy and structural stability for discrete dynamical systems”  In _Bulletin of the American Mathematical Society_ 78.6 American Mathematical Society, 1972, pp. 923–952
  * [18] N.H. Kuiper and J.W. Robbin  “Topological classification of linear endomorphisms”  In _Inventiones Mathematicae_ 19, 1973, pp. 83–106
  * [19] Clark Robinson  “Dynamical Systems: Stability, Symbolic Dynamics, and Chaos”  CRC Press, 1995
  * [20] Yuri A. Kuznetsov  “Elements of Applied Bifurcation Theory”  Springer, 2004
  * [21] R Abraham, J.E. Marsden and T Ratiu  “Manifolds, Tensor Analysis, and Applications”  Springer, 1988
  * [22] John M. Lee  “Introduction to Topological Manifolds”  Springer, 2011
  * [23] Sean Meyn and Richard L Tweedie  “Markov Chains and Stochastic Stability”  Cambridge University Press, 2009
  * [24] Panos J. Antsaklis and Anthony N. Michel  “Linear Systems”  Birkhäuser, 2006
  * [25] James D. Hamilton  “Time Series Analysis”  Princeton University Press, 1994
  * [26] Gene H. Golub and Charles F. Loan  “Matrix Computations”  John Hopkins University Press, 2013
  * [27] Alon Cohen et al.  “Online Linear Quadratic Control”  In _International Conference on Machine Learning_ , 2018, pp. 1029–1038
  * [28] Chien-Hua Lee  “Upper and lower matrix bounds of the solution for the discrete Lyapunov equation”  In _IEEE Transactions on Automatic Control_ 41.9, 1996, pp. 1338–1341