11institutetext: Institut d’Astrophysique et de Géophysique de l’Université de
Liège, Allée du 6 Août 17, 4000 Liège, Belgium 22institutetext: Observatoire
de Paris, LESIA, CNRS UMR 8109, F-92195 Meudon, France 33institutetext: Georg-
August-Universität Göttingen, Institut für Astrophysik, Friedrich-Hund-Platz
1, D-37077 Göttingen, Germany
33email<EMAIL_ADDRESS>
# Non-perturbative effect of rotation
on dipolar mixed modes in red giant stars

R-M. Ouazzani 11 2 2 M.J. Goupil 22 M-A. Dupret 11 J.P. Marques 3322

(Received , 2012; accepted , 2013)

###### Abstract

Context. The space missions CoRoT and Kepler provide high quality data that
allow us to test the transport of angular momentum in stars by the seismic
determination of the internal rotation profile.

Aims. Our aim is to test the validity of the seismic diagnostics for red giant
rotation that are based on a perturbative method and to investigate the
oscillation spectra when the validity does not hold.

Methods. We use a non-perturbative approach implemented in the ACOR code
(Ouazzani et al. 2012) that accounts for the effect of rotation on pulsations,
and solves the pulsation eigenproblem directly for dipolar oscillation modes.

Results. We find that the limit of the perturbation to first order can be
expressed in terms of the rotational splitting compared to the frequency
separation between consecutive dipolar modes. Above this limit, non-
perturbative computations are necessary but only one term in the spectral
expansion of modes is sufficient as long as the core rotation rate remains
significantly smaller than the pulsation frequencies. Each family of modes
with different azimuthal symmetry, $m$, has to be considered separately. In
particular, in case of rapid core rotation, the density of the spectrum
differs significantly from one $m$-family of modes to another, so that the
differences between the period spacings associated with each $m$-family can
constitute a promising guideline toward a proper seismic diagnostic for
rotation.

###### Key Words.:

asteroseismology - stars: interiors - stars: oscillations

##  1 Introduction

Seismic measurements of rotation profiles inside the Sun as well as stars
provide tight constraints on models of transport of angular momentum
(Pinsonneault et al., 1989; Zahn, 1992; Zahn et al., 1997; Talon & Charbonnel,
2008, and references therein). In particular, stars in late stages of
evolution, due to the highly condensed core, oscillate with non-radial modes
that have a mixed character: they behave as p modes in the envelope and as g
modes in the core. These modes, also known as mixed modes, are of particular
interest for the determination of the rotation profile throughout the star, as
they carry the signature of the star’s innermost layers and are detectable at
the surface. The CoRoT (Baglin et al., 2006) and _Kepler_ (Borucki et al.,
2010) spacecrafts have dramatically improved the quality of the available
asteroseismic data. Several recent studies reported the detection of mixed
modes that are split by rotation in a subgiant (Deheuvels et al., 2012) and in
several red giants (Beck et al., 2012; Mosser et al., 2012b, a) observed with
_Kepler_. A large number of these stars exhibit frequency spectra that show a
quite simple structure where symmetric patterns around axisymmetric modes are
easily identified. They are interpreted as multiplets of modes split by
rotation. The rotational splittings, i.e. the frequency spacing related to the
lift of degeneracy caused by rotation is then used to determine the core
rotation. The values of the corresponding splittings are quite small and the
use of the lowest order approximation to derive the splittings from stellar
models can be justified. Such studies led to the determination of unexpectedly
low central rotation frequencies (of few hundreds of nHz). These results are
in strong disagreement with the core rotation frequencies predicted by
evolutionary models, which are of the order of few dozens of $\rm\mu$Hz
(Eggenberger et al., 2012; Marques et al., 2013). They show that the transport
processes currently included in stellar models are not able to spin down the
core of red giant stars enough to explain the slowly core rotating red giants.
On the other hand, a large set of red giant stars show complex frequency
spectra (Mosser et al., 2012a), in particular with non symmetric multiplets
and therefore are likely rotating fast. Their rotation must then be
investigated with non-perturbative methods.

In this context, we first report on the relevance of using a first order
approach for the inference of rotation from seismic spectra of red giant stars
with slowly to rapidly rotating core. When not relevant, we adopt the non-
perturbative approach in order to shed light on the behaviour of splitted
mixed modes in red giants spectra.

##  2 Theoretical frequency spectra for red giants

Figure 1:  Pulsation frequencies in the inertial frame $\nu$ (in $\mu$Hz)
versus core rotation frequency $\Omega_{c}/2\pi$ (in $\mu$Hz) for dipolar
multiplets. Red curves refer to $\rm m=-1$ (prograde) modes, green curves to
$\rm m=0$ modes (axis-symmetric) , and blue curves to $\rm m=1$ (retrograde)
modes, for central rotation frequency ranging from 0 to 100 $\mu$Hz.

This work is based on the study of a model, (model M1with mass 1.3 M⊙ and
radius 3.6 R⊙), at the bottom of the red giant branch. The stellar model is
computed with the CESTAM code (Code d’Evolution Stellaire, avec Transport,
Adaptatif et Modulaire, Marques et al., 2013). Transport of angular momentum
induced by rotation is included according to Zahn (1992). The central rotation
rate in this model is ${\rm\Omega_{c}/2\pi\simeq 180\mu}$Hz, while the surface
convective region rotates at a rate of $1\,\mu$Hz. The distortion of the model
due to the centrifugal force can be neglected everywhere in red giant stars.
For the model M1, scaling laws (Kjeldsen & Bedding, 1995) give a frequency of
maximum power around $\rm\nu_{max}=289\,\rm\mu$Hz and a large separation of
$\Delta\nu=23\,\rm\mu$Hz. We then compute frequencies ranging between
$\rm\nu_{max}\pm 2\Delta\nu$. In this frequency range, the impact of the
Coriolis force remains small except in the very inner layers of the star where
it can significantly affect the modes. In order to investigate the effect of
core rotation on the frequency spectrum, we compute sets of frequencies for
model M1 for a sequence of rotation profiles. This sequence is obtained by
dividing the whole rotation profile given by CESTAM for model M1 by constant
factors. The oscillation frequencies are calculated by the non-perturbative
pulsation code ACOR (Adiabatic Code of Oscillations including Rotation,
Ouazzani et al., 2012). The eigenmodes are obtained as a result of the
coupling of spherical harmonics. In what follows, M coupling terms means
expansion on the $\ell=1,3,...,2\rm M+1$ spherical harmonics for the scalars
and the poloidal velocity component and $\ell=2,4,...,2\rm M+2$ for the
toroidal velocity component.

Figure 1 shows the frequencies of several dipolar ($\ell=1$) multiplets with
increasing rotation frequency, computed using one coupling term ($M=1$) in the
spectral expansion. Starting at low rotation rate, due to the combined action
of the Doppler effect and the Coriolis force, in the inertial frame, prograde
($m=-1$) modes are shifted towards higher frequencies, whereas retrograde
($m=+1$) modes are shifted towards lower frequencies. Crossings between modes
of different symmetry occur approximately at $\Omega_{c}/2\pi\gtrsim 8\mu$Hz
(for model M1). As a result, for core rotation frequencies above $8\mu$Hz,
modes of different $m$ are no longer gathered by original multiplets, i.e.
triplets of modes which have the same degenerate frequency without rotation.

The choice of $M=1$ coupling term in the spectral expansion is indeed
sufficient for most of the rotation profiles investigated here. This is
illustrated in Fig. 2 which shows the impact of the number of spherical
harmonics used in the eigenmodes expansion. Fig. 2 present the results of non-
perturbative calculations using three coupling terms $(M=3$). Most modes are
clearly dominated by their $\ell=1$ component (Fig. 2, top). Figure 2 (bottom)
shows that the frequency difference between calculations including one and
three spherical harmonics are smaller then the frequency difference between
two consecutive mixed modes (for $\Omega_{c}/2\pi<100\mu$Hz by two orders of
magnitude). From now on, we adopt $M=1$.

Figure 2: Top: Ratio between the contributions to kinetic energy of the
$\ell=3$ component and the $\ell=1$ component of a dominated $\ell=1$ mode.
Bottom: Difference between the pulsation frequencies computed with three
coupling terms ($M=3$) and with one coupling term ($M=1$), divided by the
frequency separation between two consecutive mixed modes. Computed for the
triplet indicated in thick lines in Fig. (1).

The trapping of modes essentially depends on their frequency (see for instance
Unno et al., 1989). As seen in Fig. 1, when $\Omega_{c}$ increases, so does
the frequency difference between the three members of a triplet and the
trapping of these members can be significantly different, to such an extent
that they end up with very different p-g nature. This is illustrated in Fig.
3, where the kinetic energy corresponding to modes of different $m$ around
$\nu_{max}$ is plotted. It gives an indication of the p-g nature of modes: the
p-dominated modes (referred as p-m modes) correspond to the minima of energy,
while the g-dominated ones (g-m modes) are associated with the maxima. The
three modes circled in black belong to the same original multiplet, and show
different p-g nature. The prograde and the retrograde modes are g-m modes,
while the $m=0$ one is a p-m mode. This change of nature induced by rotation,
which depends on the azimuthal order, can also be characterized by the number
of nodes in the p-mode cavity ($\rm n_{p}$) and in the g-mode cavity ($\rm
n_{g}$) (calculated according to the Cowling approximation, Cowling, 1941). As
shown in Fig. 3, the $\rm n_{p}$ and $\rm n_{g}$ values of the members of an
original triplet are modified differently. These gains (for $m=0,1$ modes) or
losses (for $m=-1$ modes) of nodes occur during gravito-acoustic avoided
crossings. When increasing the rotation rate, if the frequencies of two modes
of same symmetry become very close, they avoid to cross each other and
exchange nature. These avoided crossings need to be taken into account all the
more because they mainly affect the p-m modes that are the most likely to be
observed. Therefore, modes of different azimuthal order $m$ probe differently
the stellar interior. Even if they belong to the same triplet with radial
order $n$ and degree $\ell$, they can be of very different nature, and using
the rotational splitting
$\delta\omega_{n,\ell}=(\omega_{n,\ell,m}-\omega_{n,\ell,-m})/2m$ to determine
the rotation rate is therefore questionable for red giants.

##  3 Slow to moderate core rotating red giants

For models with $\Omega_{c}/2\pi\leq$ 8 $\rm\mu$Hz (see Fig. 1), the
frequencies behave linearly with respect to $\Omega_{c}$. In this range, the
sectorial modes ($\ell=\,\mid m\mid$) are symmetrically distributed around the
axisymmetric modes, and have the same radial order.

Figure 3:  For a central rotation frequency of 20$\rm\mu$Hz, kinetic energy of
modes around a p-m mode as a function of the mode frequency. The three modes
circled in black belong to the same original triplet. Inside the parentheses
are indicated their number of nodes in the p cavity and in the g cavity.

In this case, a first order perturbation approach approximates well the
effects of rotation on the mode frequencies. The pulsation $\omega_{n,\ell,m}$
is then given by
$\omega_{n,\ell,m}\,=\,\omega_{n,\ell,0}\,+\,m\delta\omega_{n,\ell}^{1}$
(Cowling & Newing 1949 and Ledoux 1949) where the first order rotational
splitting $\delta\omega_{n,\ell}^{1}$ is expressed as a weighted measure of
the star rotation rate:
$\delta\omega_{n,\ell}^{1}\,=\,\int_{0}^{R}\,K_{n,\ell}(r)\,\Omega(r)\,\rm
dr$, where $K_{n,\ell}(r)$ are the rotational kernels of the modes. They
depend on the equilibrium structure and on the eigenfunctions of the
unperturbed modes (see e.g. Goupil, 2009, and references therein). This
formulation relies on two main assumptions. First, only the Coriolis force is
significant and it is enough to account for its contribution to order
$O(\Omega/\omega_{n,\ell})$. Second, it is derived from a variational
principle (see e.g. Lynden-Bell & Ostriker, 1967), which requires that the
eigenfunction of the mode perturbed by rotation is close to the eigenfunction
of the corresponding unperturbed mode. This is the case for model M1 for
$\Omega_{c}/2\pi<8\leavevmode\nobreak\ \mu$Hz.

Between 8 and 20 $\rm\mu$Hz (Fig. 1), the mode frequencies seem to behave
linearly with $\Omega_{c}$, but are no longer gathered by original multiplet.
Figure 4 (top) displays the apparent rotational splittings, taking half the
difference between the closest modes of opposite azimuthal order $m=\pm 1$.
The values taken by the apparent rotational splittings range from few hundreds
nHz to 2 $\rm\mu$Hz. This is of the order of the splittings measured in
observed spectra. Because the apparent multiplets occur by ‘accident’, the
curve in Fig. 4 (top) does not –and usually cannot– follow the particular ‘V’
pattern around the p-m modes as observed for some slowly rotating red giants
(see for observations Beck et al. 2012, Mosser et al. 2012b). This ‘V’ pattern
has the same origin as the modulation seen for mode inertia (Goupil et al.
2013) that is due to the trapping of modes (see Dziembowski et al., 2001;
Dupret et al., 2009). The ‘V’ pattern can then constitute a strong indication
that the selection of multiplet has been done correctly in observed spectra.
Figure 3 confirms that modes which belong to the same original multiplet
(circled in black) are found in different apparent multiplets, and show quite
different p-g nature. Fig 4 (bottom) compares the original splittings –i.e.
half the differences between frequencies of modes from the same original
triplet with opposite $m$– with the splittings given by the first order
approach. The g-m modes (e.g. at 277 $\mu$Hz and 303 $\mu$Hz) have the largest
splittings. For such modes, both computations give very close results (with a
discrepancy of $\sim 4\%$). The p-m modes (e.g. at 290 $\mu$Hz and 315
$\mu$Hz) have the smallest splittings and show the expected ‘V’ shape
variation. For these modes, the discrepancies can reach $30\%$. This is due to
the trapping of these modes which differs from one $m$ to the other. Note that
original splittings including three terms ($M=3$) in the non-perturbative
calculations give relative differences lower than $10^{-6}$ when compared with
results of computations using two $Y_{\ell}^{m}$ (Fig. 2).

Figure 4:  For the model M1 and a central rotation frequency of $20\rm\mu$Hz,
Top: spacing between two closest modes of opposite azimuthal order $\rm m=\pm
1$, i.e. apparent rotational splitting. Bottom: spacing between prograde and
the retrograde modes which belong to the same original multiplet (grey,
crosses) and splittings computed with the first order perturbative method
(green, plus signs).

Therefore, for a rotation between 8 and 20 $\rm\mu$Hz, the fact that the
multiplets are correctly selected can be assumed if the splitting given by
apparent triplets do follow the ‘V’ pattern. If so, the first order approach
provides an order of magnitude estimate of the core rotation rate, but non-
perturbative modeling is required for accurate quantitative conclusions on the
rotational profile from p-m modes.

##  4 Rapid core rotating red giants

For $\Omega>20\mu$Hz, the mode frequencies no longer behave linearly with the
rotation rate (see Fig. 1). This is due not only to higher order effects that
come into play, but also to the trapping of modes that is modified by
rotation. Modes which belong to the same original multiplet do not have the
same radial distribution. Under these circumstances, the notion of rotational
splitting as defined by $\delta\omega_{n,\ell}^{1}$ is no longer relevant, and
cannot be simply related to the rotation profile. Figure 5 (top) shows the
kinetic energy of pulsating modes for a core rotation rate of 140 $\rm\mu$Hz.
The large separation, as measured by the difference between consecutive p-m
modes with minimum kinetic energies, is conserved. The higher kinetic energy
of prograde modes indicates that they are more of g nature than retrograde or
axisymmetric ones. This is due to the shift in frequency induced by rotation
that brings, in the same frequency range, modes which belong to very different
parts of the zero rotation spectrum (see Fig. 1).

Figure 5 (top) shows that, in the same range of frequency, there are more
$m=-1$ modes than $m=0$ modes, and in turn, there are more $m=0$ modes than
$m=+1$ ones. In order to highlight this difference in distribution, the period
spacings are plotted for each $m$ value in Fig. 5 (bottom). In this diagram,
the families of modes of different azimuthal orders clearly show different
values of period spacing, ranging from the lowest for the prograde modes to
the highest for the retrograde ones. Note that this phenomenon appears only
for a rotation high enough that the distributions of the three $m$-families of
modes become clearly different ($\Omega_{c}\sim 100\mu$Hz for model M1). If
one is able to measure three different values of period spacings in a observed
spectrum, one is able to identify the values of the azimuthal order. Moreover,
based on our study of synthetic spectra, we find that the differences between
the values of period spacings increase with $\Omega_{c}$, and therefore can
constitute a guideline toward a constraint on the rotation profile.

##  5 Discussions and conclusions

Figure 5: For a central rotation frequency of 140 $\rm\mu$Hz, top: pulsation
modes’ kinetic energy with respect to their frequency. Bottom: period spacing
with respect to the pulsation frequency.

Some observed spectra of red giant stars show structures with nearly symmetric
spacings around axisymmetric modes. For this kind of spectra, we found that
one can figure out if these apparent multiplets correspond to the original
splittings by plotting the apparent splittings as a function of the
axisymmetric modes frequencies. If the resulting curve follows the ‘V’
pattern, then it is reasonnable to assume that the triplets are correctly
selected, and quantitative information on rotation can be derived from the
observed splittings. This corresponds to a low rotation regime when the
splitting is smaller than the frequency spacing between two consecutive modes
of same $m$, i.e. for ${\rm\Omega_{c}/(2\pi)/(\Delta P/P^{2})}\lesssim 2$.

If the apparent splittings do not follow the ‘V’ pattern, the triplets
overlap. The star is therefore located in the regime of moderate rotation. The
frequency differences between the original triplet members are large enough
that the trapping properties differ from one member to another. However, if it
is still possible to select correctly the original triplets (as done in Mosser
et al., 2012a), the first order approach could be used for a first guess of
the core rotation rate. Non-perturbative calculations are nevertheless
required for a precise determination, particularly because of the avoided
crossings that p-m modes undergo. This rotation regime corresponds to values
${\rm\Omega_{c}/(2\pi)/(\Delta P/P^{2})}$ between 2 and 5.

Finally, some observed spectra of red giants do not show any regular or close
to regular structures. We expect these spectra to correspond to rapid
rotators. In this case, modes of different azimuthal orders have a very
different nature. The concept of rotational splitting of modes with the same
radial distribution and different $m$ is no longer relevant, and one should
consider separately sub-spectra associated to each value of $m$. Provided that
the rotation is high enough to give rise to clearly different distributions
with respect to $m$, the differences between the period spacings associated
with $m$ allows to identify the azimuthal order, and thereby offers the
promising opportunity of deriving a proper seismic diagnostic. That
corresponds to the very rapid rotation case i.e. for
${\rm\Omega_{c}/(2\pi)/(\Delta P/P^{2})\gtrsim 20}$. An intermediate case
remains, where the rotational splitting is no longer a relevant seismic
diagnostic, but where the rotation is too slow to allow to distinguish three
different period spacings associated to the $m$-families of modes (for
${\rm\Omega_{c}/(2\pi)/(\Delta P/P^{2})}$ between 5 and 20). Establishing a
diagnostic in this regime constitutes an important issue which needs to be
adressed in a forthcoming paper.

It can be surprising that when $\Omega/(2\pi)/(\Delta P/P^{2})<20$, on the one
hand the perturbative approach gives inaccurate results, but on the other hand
non-perturbative models where the spectral expansions are limited to one term
(plus the toroidal component) give accurate enough results. This can be
understood with the help of the following equation:
$\Omega(r)/2\pi\sim\Delta\nu\ll\nu_{max}$. In the frequency range of interest
here, the rotation frequency remains significantly lower than $\nu_{max}$
throughout the whole star, and thus lower than the pulsation frequencies.
Hence, the Coriolis terms are smaller than other terms in the equation of
motion and the coupling of $\ell=1$ modes with $\ell=3,5,...$ components is
small. However the rotational splitting is not small compared to the large
separation $\Delta\nu$. This implies that the trapping of each member of a
same triplet can be very different, as illustrated in Fig.3. Using the
variational principle to model rotational splittings is thus not justified
because the eigenfunctions are too different. The only way to properly model
the effect of rotation on oscillations is to solve separately the differential
equations associated to each member of the multiplet.

To sum up, the major impact of moderate rotation on red giants’ spectra is the
modification of the trapping that depends on $m$. Let us precise that here the
cavities are not considered to be modified by rotation, but it is the way
modes probe them which differs from one member of a multiplet to another. In
such a case, the first order perturbative approach gives inaccurate results
and the $m$-families of modes carry different information on the stellar
interior. Only methods which compute these $m$-sub-spectra independently, i.e.
non-perturbative methods, are appropriate for studing the moderate to rapidly
core rotating red giants.

###### Acknowledgements.

RMO thanks Benoît Mosser for fruitful discussions. The authors thank the
referee for his comments that helped to improve the manuscript.

## References

  * Baglin et al. (2006) Baglin, A., Auvergne, M., Barge, P., et al. 2006, in ESA Special Publication, Vol. 1306, ESA Special Publication, ed. M. Fridlund, A. Baglin, J. Lochard, & L. Conroy, 33
  * Beck et al. (2012) Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Cowling (1941) Cowling, T. G. 1941, MNRAS, 101, 367
  * Cowling & Newing (1949) Cowling, T. G. & Newing, R. A. 1949, ApJ, 109, 149
  * Deheuvels et al. (2012) Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19
  * Dupret et al. (2009) Dupret, M.-A., Belkacem, K., Samadi, R., et al. 2009, A&A, 506, 57
  * Dziembowski et al. (2001) Dziembowski, W. A., Gough, D. O., Houdek, G., & Sienkiewicz, R. 2001, MNRAS, 328, 601
  * Eggenberger et al. (2012) Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4
  * Goupil (2009) Goupil, M. J. 2009, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 765, The Rotation of Sun and Stars, 45–99
  * Goupil et al. (2013) Goupil, M. J., Mosser, B., Marques, J. P., et al. 2013, A&A, 549, A75
  * Kjeldsen & Bedding (1995) Kjeldsen, H. & Bedding, T. R. 1995, A&A, 293, 87
  * Ledoux (1949) Ledoux, P. 1949, Memoires of the Societe Royale des Sciences de Liege, 9, 3
  * Lynden-Bell & Ostriker (1967) Lynden-Bell, D. & Ostriker, J. P. 1967, MNRAS, 136, 293
  * Marques et al. (2013) Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74
  * Mosser et al. (2012a) Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012a, A&A, 548, A10
  * Mosser et al. (2012b) Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012b, A&A, 540, A143
  * Ouazzani et al. (2012) Ouazzani, R.-M., Dupret, M.-A., & Reese, D. R. 2012, A&A, 547, A75
  * Pinsonneault et al. (1989) Pinsonneault, M. H., Kawaler, S. D., Sofia, S., & Demarque, P. 1989, ApJ, 338, 424
  * Talon & Charbonnel (2008) Talon, S. & Charbonnel, C. 2008, A&A, 482, 597
  * Unno et al. (1989) Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars
  * Zahn (1992) Zahn, J. 1992, A&A, 265, 115
  * Zahn et al. (1997) Zahn, J.-P., Talon, S., & Matias, J. 1997, A&A, 322, 320


Using $M^+$ to denote the Moore-Penrose Inverse of a matrix $M$, we can also write $\bm{u}^*_{c} :=  P^{+}_H A^T e_j = P^\top_H W^{-1}_H A^T e_j  $, where $W_H = P_H P_H^\top$.

Ignoring the constant term $(C^TC)_{jj}$, we can reparametrize $\bm{u}= -\bm{u}^*_c + \bm{u}'$ where $\bm{u}'\in\mathrm{null}{N}(P_H)$ and describe (<ref>) as the optimization problem:

\begin{align}\label{eq:sj}
S_j \quad &:= \min\limits_{\bm{u}' \in \mathrm{null}(P_H(A,B))} \Big\|\begin{bmatrix}\bm{C}& 0\\ 0 & \bm{D}\end{bmatrix}\begin{bmatrix}\bm{G}_u(A,B) \\ I \end{bmatrix}(\bm{u}'-\bm{u}_c^*(A,B))\Big\|^2_2
\end{align}

The least-squares problem.

Let $H$ denote the FIR-Horizon of the problem, then define the matrices

\begin{align}\label{eq:gab}
&\bm{G}_w(A) = \begin{bmatrix}I &0 & 0 & \dots & 0\\
A &I & 0 & \dots & 0\\
A^2 &A & I & \dots & 0\\
& \dots & \dots & & \\
A^{H-1} & A^{H-2} & A^{H-3} & \dots & I\\
\end{bmatrix}, &\bm{G}_u(A,B) = \begin{bmatrix}B &0 & 0 & \dots & 0\\
AB &B & 0 & \dots & 0\\
A^2B &AB & B & \dots & 0\\
& \dots & \dots & & \\
A^{H-1}B & A^{H-2}B & A^{H-3}B & \dots & B\\
\end{bmatrix}.
\end{align}

and denote $P_i(A,B)$ as the $i$th block matrix row of $\bm{G}_u(A,B)$:

\begin{align}\label{eq:pdef}
P_i(A,B) = [A^{i-1}B, A^{i-2}B, \dots, B, 0, \dots,0]
\end{align}

$\bm{G}_u(A,B)$ can be written as $\bm{G}_u(A,B) = \bm{G}_w(A)(I_H \otimes B)$, where $I_H$ is the identity matrix in $\R^H$. Let $Z \in \mathbb{R}^{H \times H}$ be defined as the nilpotent matrix

\begin{align}
Z = \begin{bmatrix} \bm{0}_{H-1 \times 1}  & I_{H-1} \\
0& \bm{0}_{1 \times H-1} \end{bmatrix},
\end{align}

and notice it's psuedo-inverse is $Z^+ = Z^\top$. Using $Z$, it is easy to verify that $\bm{G}_w(A)$ can be expressed as:

\begin{align}
\bm{G}_w(A) = \left(I_H - Z^+ \otimes A\right)^{-1}
\end{align}

Recall our derivation of the reduced problem $S_j$ (<ref>) for a fixed $j$:

\begin{align}\label{eq:sj}
S_j \quad &:= \min\limits_{\bm{u} \in \mathrm{null}(P_H(A,B))} \Big\|\begin{bmatrix}\bm{C}& 0\\ 0 & \bm{D}\end{bmatrix}\begin{bmatrix}\bm{G}_u(A,B) \\ I \end{bmatrix}(\bm{u}-\bm{u}_c(A,B))\Big\|^2_2
\end{align}

Where $\bm{u}_{c}$ denotes $\bm{u}_{c} :=  P^{+}_H A^H e_j = P^\top_HW^{-1}_HA^H e_j$. Let $\bm{u}^*(A,B)$ be a minimizer of the above problem for fixed $A,B$, we are interested in the SLS solutions
$$\phi^{*j}(A,B) := \begin{bmatrix}\phi^{*j}_x(A,B) \\ \phi^{*j}_u(A,B) \end{bmatrix} = \begin{bmatrix}\bm{G}_u(A,B)\\I \end{bmatrix}(\bm{u}^*(A,B)-\bm{u}_c(A,B))$$

and how these solutions are perturbed with changes in $A,B$.

For the rest of the discussion, we will drop mentioning the explicit dependence on $(A,B)$ to reduce the notational burden.
First, we (over-)parametrize $\bm{u}$ as
$$\bm{u} = (I-P^+_H P_H)\bm{\eta}, $$
to cast the above problem into an unconstrained one:

\begin{align}\label{eq:sj}
S_j \quad &:= \min\limits_{\bm{\eta}} \Big\|\underbrace{\begin{bmatrix}\bm{C}& 0\\ 0 & \bm{D}\end{bmatrix}\begin{bmatrix}\bm{G}_u(A,B) \\ I \end{bmatrix}(I-P^+_HP_H)}_{\F}\bm{\eta}- \underbrace{\begin{bmatrix}\bm{C}& 0\\ 0 & \bm{D}\end{bmatrix}\begin{bmatrix}\bm{G}_u(A,B) \\ I \end{bmatrix}\bm{u}_c(A,B)}_{\g}\Big\|^2_2
\end{align}

The unique min-norm solution $\eta^*$ to the above problem is $\eta^* = \F^+\g$ and therefore the optimal solution $\phi^*$ takes the form

\begin{align}
\phi^* = \begin{bmatrix}\bm{C}^{-1}& 0\\ 0 & \bm{D}^{-1}\end{bmatrix}(\F \F^+\g - \g) = \begin{bmatrix}\bm{C}^{-1}& 0\\ 0 & \bm{D}^{-1}\end{bmatrix}\underbrace{(\F \F^+ - I)\g}_{\nu^*} =: \begin{bmatrix}\bm{C}^{-1}& 0\\ 0 & \bm{D}^{-1}\end{bmatrix} \nu^*
\end{align}

§.§ Local lipshitzness of $\cl{H}_2$-optimal closed-loop operators

Here, we perform perturbation analysis on the term $\nu^* = (\F \F^+ - I)\g$. Throughout the discussion, we will make frequent use of the following identities:

For arbitrary matrices $X,Y\in\R^{n \times m}$ and $A, B \in \R^{n \times n}$  holds:

* $A^k_1-A^k_2 = \sum^{k-1}_{j=0} A^{k-1-j}_1(A_1-A_2)A^j_2$

* $XX^+-YY^+ = (I-XX^+)(X-Y)Y^+ + \left[(I-YY^+)(X-Y)X^+\right]^\top$

* If $A$ and $B$ are invertible, then $A^{-1}-B^{-1}= A^{-1}(B-A)B^{-1}$.

The following is a corollary from Theorem 4.1 in [Wedin, 1973]:

Let $X$ and $Y$ be matrices with equal rank, let $\|\:\cdot\:\|_2$ denote the induced 2-norm and $\|\:\cdot\:\|_F$ denote the Frobenius norm. The following inequalities hold:

\begin{align*}
\|X^+-Y^+\|_2 &\leq \varphi \|X^+\|_2\|Y^+\|_2 \|X-Y\|_2\\
\|X^+-Y^+\|_F &\leq \sqrt{2} \|X^+\|_2\|Y^+\|_2 \|X-Y\|_F
\end{align*}

where $\varphi = \tfrac{1 + \sqrt{5}}{2}$ denotes the golden ratio constant.

Next we present the core theorem of the perturbation analysis: Given two arbitrary controllable systems $(A_1,B_1)$ and $(A_2,B_2)$, thm:mainsensitivity bounds the worst-case difference in solutions $\|\phi^{*j}_1-\phi^{*j}_2\|_2$ in terms of the differences in parameters space $\|A_1-A_2\|_{F}$ and $\|B_1-B_2\|_2$ between both systems. This result is the first global (considering arbitrary pairs of $A_1,A_2$ and $B_1,B_2$) perturbation bound  for $\cl{H}_2$-optimal control with SLS.

Let $C,D \succ 0$, let $(A_1,B_1)$ and $(A_2,B_2)$ be two controllable pairs of system matrices with FIR horizon $H$ and let $\phi^{*j}_1$ and $\phi^{*j}_2$ be the corresponding SLS-solutions to the subproblem $S_j$. Then, it holds that:

\begin{align}
\|\phi^{j*}_1-\phi^{j*}_2\|_2 \leq \Gamma_A\|A_1-A_2\|_F + \Gamma_B \|B_1-B_2\|_F
\end{align}

where the Lipshitz-constants $\Gamma_A,\Gamma_B$ stand for

\begin{align*}
\Gamma_A &= \kappa_{CD}\Gamma'_1 + \kappa_{CD}\Gamma'_2 \|B_1\|_2\|\bm{G}_w(A_1)\|_2, &&\kappa_{CD} = \frac{\max\{\sigma_{max}(C),\sigma_{max}(D)\}}{\min\{\sigma_{min}(C),\sigma_{min}(D)\}}  \\
\Gamma_B &= \kappa_{CD}\Gamma'_2 \|\bm{G}_w(A_2)\|_2 &&
\end{align*}

and $\Gamma'_1$ and $\Gamma'_2$ are defined as:

\begin{align*}
\Gamma'_1 &= \alpha_{H,1} \alpha_{H,2}H(1+\left\|\bm{G}_{u,2}\right\|_2)\|P^+_{H,2}\|_2  \\
\Gamma'_2 &= \alpha_{H,1} \|P^+_{H,1}\|_2 \left(1 + \bm{\varphi} \|P^+_{H,2}\|_2 +\bm{\varphi} \|P^+_{H,2}\|_2\left\|\bm{G}_{u,2}\right\|_2 \right) + \|\g\|_2(\|\F^+_1\|_2 + \|\F^+_2\|_2) + \dots \\
\notag & \quad + \bm{\varphi}\|\g\|_2(\|\F^+_1\|_2 + \|\F^+_2\|_2)\|P^+_{H,1}\|_2 \|P^+_{H,2}\|_2(\|P_{H,1}\|_2 + \|P_{H,2}\|_2)(1+\|\bm{G}_{u,1}\|_2).
\end{align*}

and $\varphi=\tfrac{1+\sqrt{5}}{2}$ is the golden ratio.

Recall the identities of lem:diffs. Write $\nu^*_1-\nu^*_2$ as

\begin{align}
\notag    \nu^*_1-\nu^*_2 &= (\F_1\F^+_1-I)(\g_1-\g_2) + (\F_1\F^+_1 - \F_2\F^+_2)\g_2 \\
\label{eq:nu1bound}   \|\nu^*_1-\nu^*_2\|_2&\leq \|\g_1-\g_2\|_2 + \|\F_1\F^+_1 - \F_2\F^+_2\|_2\|\g\|_2,
\end{align}

where used the fact that $(\F_1\F^+_1-I)$ is a projection and therefore $\|\F_1\F^+_1-I\|_2 =1$. Rewrite $\F_1\F^+_1 - \F_2\F^+_2$ as

$$ (I-\F_1\F^+_1)(\F_1-\F_2)\F^+_2 + \left[(I-\F_2\F^+_2)(\F_1-\F_2)\F^+_1\right]^\top $$
to conclude that

\begin{align}
\|\F_1\F^+_1 - \F_2\F^+_2\|_2 \leq \|\F_1-\F_2\|_2(\|\F^+_1\|_2 + \|\F^+_2\|_2).
\end{align}

Substitution into (<ref>) yields:

\begin{align}
\label{eq:nu2bound}   \|\nu^*_1-\nu^*_2\|_2 &\leq \|\g_1-\g_2\|_2 + \|\F_1-\F_2\|_2(\|\F^+_1\|_2 + \|\F^+_2\|_2)\|\g\|_2,
\end{align}

* Bounding $\|\F_1-\F_2\|_2$: Rewrite $\F_1-\F_2$ as

\begin{align}
\begin{bmatrix}\bm{C}^{-1}& 0\\ 0 & \bm{D}^{-1}\end{bmatrix}(\F_1-\F_2) &=  \begin{bmatrix}\bm{G}_{u,1} \\ I \end{bmatrix}(I-P^+_{H,1}P_{H,1}) -\begin{bmatrix}\bm{G}_{u,2} \\ I \end{bmatrix}(I-P^+_{H,2}P_{H,2}) \\
&=  \begin{bmatrix}\bm{G}_{u,1} \\ I \end{bmatrix}(P^+_{H,2}P_{H,2}-P^+_{H,1}P_{H,1}) + \begin{bmatrix}\bm{G}_{u,1} - \bm{G}_{u,2} \\ 0 \end{bmatrix}(I-P^+_{H,2}P_{H,2})
\end{align}

From the above we can derive the inequality:

\begin{align}\label{eq:F12bound}
\frac{\|\F_1-\F_2\|_2}{\max\{\|C\|_2,\|D\|_2\}} &\leq (1+\|\bm{G}_{u,1}\|_2)\|P^+_{H,2}-P^+_{H,1}\|_2(\|P_{H,1}\|_2 + \|P_{H,2}\|_2) + \|\bm{G}_{u,1}-\bm{G}_{u,2}\|_2
\end{align}

Now we will use the result thm:wedin to bound $\|P^+_{H,2}-P^+_{H,1}\|_2$ as

\begin{align}
\|P^+_{H,2}-P^+_{H,1}\|_2 \leq \bm{\varphi} \|P^+_{H,1}\|_2 \|P^+_{H,2}\|_2 \|P_{H,2}-P_{H,1}\|_2
\end{align}

Furthermore, noticing $P_{H,2}-P_{H,1} = [\bm{0},\dots,\bm{0},\bm{I}_n](\bm{G}_{u,2}-\bm{G}_{u,1})$ we can conclude

\begin{align}
\|P^+_{H,2}-P^+_{H,1}\|_2 \leq \bm{\varphi} \|P^+_{H,1}\|_2 \|P^+_{H,2}\|_2 \|\bm{G}_{u,2}-\bm{G}_{u,1}\|_2.
\end{align}

We combine this into (<ref>) to obtain

\begin{align}\label{eq:F12bound-2}
\notag &\frac{\|\F_1-\F_2\|_2}{\max\{\|C\|_2,\|D\|_2\}}\\
\leq &\left( 1+ \bm{\varphi} \|P^+_{H,1}\|_2 \|P^+_{H,2}\|_2 (1+\|\bm{G}_{u,1}\|_2)(\|P_{H,1}\|_2 + \|P_{H,2}\|_2) \right)\|\bm{G}_{u,1}-\bm{G}_{u,2}\|_2
\end{align}

* Bounding $\|\g_1-\g_2\|_2$:
Introduce the constant $\alpha_H := \max_{0\leq k\leq H} \|A^k\|_2$ and observe that $\|A^H_1-A^H_2\|_2$ can be bounded as:

\begin{align}
\|A^H_1-A^H_2\|_2 = \|\sum^{H-1}_{j=0} A^{H-1-j}_1(A_1-A_2)A^j_2\| \leq H \alpha_{H,1} \alpha_{H,2} \|A_1-A_2\|_2
\end{align}

We can rewrite $\g_1-\g_2$ as

\begin{align}
\begin{bmatrix}\bm{C}^{-1}& 0\\ 0 & \bm{D}^{-1}\end{bmatrix}(\g_1-\g_2) &= \begin{bmatrix}\bm{G}_{u,1} \\ I \end{bmatrix}P^+_{H,1}A^H_1e_j - \begin{bmatrix}\bm{G}_{u,2} \\ I \end{bmatrix}P^+_{H,2}A^H_2e_j \\
&=\begin{bmatrix}(\bm{G}_{u,1}-\bm{G}_{u,2}) \\ 0 \end{bmatrix}P^+_{H,1}A^H_1e_j + \begin{bmatrix}\bm{G}_{u,2} \\ I\end{bmatrix}(P^+_{H,1}-P^+_{H,2})A^H_1e_j \\
\notag&\quad \dots+ \begin{bmatrix}\bm{G}_{u,2} \\ I \end{bmatrix}P^+_{H,2}(A^H_1-A^H_2)e_j
\end{align}

and obtain the bound:

\begin{align}
\frac{\|\g_1-\g_2\|_2}{\max\{\|C\|_2,\|D\|_2\}} &\leq \alpha_{H,1}\left\| \bm{G}_{u,1}-\bm{G}_{u,2}\right\|_2 \|P^+_{H,1}\|_2+ \alpha_{H,1}(1+\left\|\bm{G}_{u,2}\right\|_2)\left\|P^+_{H,1}-P^+_{H,2}\right\|_2 \\
&\quad \dots + \alpha_{H,1} \alpha_{H,2}H(1+\left\|\bm{G}_{u,2}\right\|_2)\|P^+_{H,2}\|_2\|A_1-A_2\|_2 \\
&\leq  \alpha_{H,1} \|P^+_{H,1}\|_2 \left(1 + \bm{\varphi} \|P^+_{H,2}\|_2 +\bm{\varphi} \|P^+_{H,2}\|_2\left\|\bm{G}_{u,2}\right\|_2 \right)\left\| \bm{G}_{u,1}-\bm{G}_{u,2}\right\|_2 \\
&\quad \dots + \alpha_{H,1} \alpha_{H,2}H(1+\left\|\bm{G}_{u,2}\right\|_2)\|P^+_{H,2}\|_2\|A_1-A_2\|_2
\end{align}

We get the bound

\begin{align}
\frac{\|\nu^*_1-\nu^*_2\|_2}{\max\{\|C\|_2,\|D\|_2\}} \leq \Gamma'_1 \|A_1-A_2\|_2 + \Gamma'_2 \|\bm{G}_{u,1}-\bm{G}_{u,2}\|_2
\end{align}

where $\Gamma'_1$ and $\Gamma'_2$ are the constants:

\begin{align}
\Gamma'_1 &= \alpha_{H,1} \alpha_{H,2}H(1+\left\|\bm{G}_{u,2}\right\|_2)\|P^+_{H,2}\|_2  \\
\Gamma'_2 &= \alpha_{H,1} \|P^+_{H,1}\|_2 \left(1 + \bm{\varphi} \|P^+_{H,2}\|_2 +\bm{\varphi} \|P^+_{H,2}\|_2\left\|\bm{G}_{u,2}\right\|_2 \right) + \|\g\|_2(\|\F^+_1\|_2 + \|\F^+_2\|_2) + \dots \\
\notag & \quad + \bm{\varphi}\|\g\|_2(\|\F^+_1\|_2 + \|\F^+_2\|_2)\|P^+_{H,1}\|_2 \|P^+_{H,2}\|_2(\|P_{H,1}\|_2 + \|P_{H,2}\|_2)(1+\|\bm{G}_{u,1}\|_2)
\end{align}

Using lem:Gwu, we obtain the final bound:

\begin{align}\label{eq:Gamma12}
\|\phi^*_1 - \phi^*_2\|_2 \leq  \kappa_{CD}\|\nu^*_1-\nu^*_2\|_2 \leq \Gamma_A \|A_1-A_2\|_2 + \Gamma_B \|B_1-B_2\|_2
\end{align}

with the constants $\Gamma_A,\Gamma_B$ defined as:

\begin{align}\label{eq:GammaAB}
\Gamma_A &= \kappa_{CD}\Gamma'_1 + \kappa_{CD}\Gamma'_2 \|B_1\|_2\|\bm{G}_w(A_1)\|_2 \|\bm{G}_w(A_2)\|_2 \\
\Gamma_B &= \kappa_{CD}\Gamma'_2 \|\bm{G}_w(A_2)\|_2
\end{align}

§.§ Global lipshitzness of $\cl{H}_2$-optimal closed-loop operators over compact sets $\cl{S}$

For the next sections we will assume that we are given some fixed compact set $\cl{S}$ of controllable systems. Recall the basic implications of lem:baseCtrl; we define accordingly the FIR horizon $H$ and the constants $\overline{\sigma}^w$, $\underline{\sigma}^w$, $\overline{\sigma}^u$, $\underline{\sigma}^u$ and assume that these quantities are all known.

This section derives a global Lipshitz bound for $\cl{H}_2$-optimal SLS solutions over a compact set of controllable systems $\cl{S}$. As a starting point we consider the previous theorem thm:mainsensitivity. Our main proof strategy is to derive global bounds on the constants $\Gamma_A$ and $\Gamma_B$. We proceed with a collection lemmas bounding individual terms in the equations (<ref>), (<ref>).

§.§.§ Auxiliary Lemmas

For any pair of system matrices $(A_1,B_1)$ and $(A_2,B_2)$ (with compatible dimensions) holds

\begin{align}
\|\bm{G}_w(A_1)-\bm{G}_w(A_2)\|_2 & \leq \|\bm{G}_w(A_1)\|_2 \|\bm{G}_w(A_2)\|_2 \| A_1 -  A_2\|_2 \\
\notag \|\bm{G}_u(A_1,B_1)-\bm{G}_u(A_2,B_2)\|_2& \leq \|B_1\|_2\|\bm{G}_w(A_1)\|_2 \|\bm{G}_w(A_2)\|_2 \| A_1 -  A_2\|_2 + \|\bm{G}_w(A_2)\|_2\|B_1- B_2\|_2
\end{align}

Using lem:diffs we can write
Using $\bm{G}_u(A,B) = \bm{G}_w(A)(I_H \otimes B)$ and lem:diffs we can write $\bm{G}_{u,1}-\bm{G}_{u,2}$ as

\begin{align}
\bm{G}_{u,1}-\bm{G}_{u,2} &= \bm{G}_w(A_1)(I_H \otimes B_1) - \bm{G}_w(A_2)(I_H \otimes B_2) \\
&= \left(\bm{G}_w(A_1)-\bm{G}_w(A_2)\right)(I_H \otimes B_1) + \bm{G}_w(A_2)\left( I_H \otimes (B_1- B_2)\right)
\end{align}

It holds that

\begin{align}
\bm{G}_w(A_1)-\bm{G}_w(A_2) &= \bm{G}_w(A_1)(\bm{G}_w(A_2)^{-1}-\bm{G}_w(A_1)^{-1})\bm{G}_w(A_2) \\
&= \bm{G}_w(A_1)(Z^+ \otimes (A_1 -  A_2))\bm{G}_w(A_2)
\end{align}

which leads to the bound

\begin{align}
\left\| \bm{G}_w(A_1)-\bm{G}_w(A_2) \right\|_2 \leq \|\bm{G}_w(A_1)\|_2 \|A_1-A_2\|_2 \|\bm{G}_w(A_2)\|_2
\end{align}

In total, we need to global bounds on the quantities $\|\bm{G}_{u}\|_2$,$\|\bm{G}_{w}\|_2$, $\|P^+_{H}\|_2$, $\|P_{H}\|_2$, $\|\F^+\|_2$, $\|\g\|_2$.

Let $(A,B)$ be pair of fixed system matrices, let $\bm{G}_{u}(A,B)$, $\bm{G}_w(A)$ be the matrices defined in (<ref>), and let $W^u_H = \sum^{H-1}_{i=0} A^iBB^\top A^{i\top}$, $W^w_H = \sum^{H-1}_{i=0} A^i A^{i\top}$  be the $H$th controllability grammian w.r.t to the input $u$ and the distrubance $w$, respectively. Then it holds:

\begin{align}
&\|\bm{G}_{u}(A,B)\|_2 \leq \sqrt{H \sigma_{max}(W^u_H(A,B))} && \|\bm{G}_{w}(A)\|_2 \leq \sqrt{H \sigma_{max}(W^w_H(A))}
\end{align}

$\|\bm{G}_{u}\|_2$ is defined as $\|\bm{G}_{u}\|^2_2 := \max\limits_{\|u\|_2 = 1} \|\bm{G}_u \bm{u}\|^2_2$, by decomposing $\bm{u} = [u^\top_0,\dots,u^\top_{H-1}]^\top$ we can rewrite this as

\begin{align}
\|\bm{G}_{u}\|^2_2 &= \max\limits_{\|u\|_2 = 1} \left\| \begin{bmatrix} B u_0 \\ ABu_0 + Bu_1\\ \dots \\ A^{H-1}Bu_0 + \dots + Bu_{H-1} \end{bmatrix} \right\|^2_2 = \max\limits_{\|u\|_2 = 1} \sum^{H}_{k=1} \|P_k \bm{u}\|^2_2 \\
&\leq  \sum^{H}_{k=1} \max\limits_{\|u\|_2 = 1} \|P_k \bm{u}\|^2_2 =\sum^{H}_{k=1}  \|P_k\|^2_2 \leq H \|P_H\|^2_2 \leq H \|W^u_H\|_2
\end{align}

Where we used the fact that $\|P_k\|^2_2$ increases in $k$ and that $\|P_k\|^2_2$ is equal to the induced $2$-norm of the corresponding controllabillity grammian $W^u_k = \sum^{k-1}_{i=0} A^iBB^\top A^{i\top}$. Thus, we obtain the bound

$$\|\bm{G}_{u}(A,B)\|_2 \leq \sqrt{H \sigma_{max}(W^u_H(A,B))}, $$

and the bound on $\|\bm{G}_{w}(A)\|_2$ follows in the same way.

Let $(A,B)$ be pair of $H$-controllable fixed system matrices, let $P_H(A,B)$ be the matrix defined in (<ref>), and let $W^u_H = \sum^{H-1}_{i=0} A^iBB^\top A^{i\top}$ be the $H$th controllability grammian w.r.t to the input $u$. Then, the induced $2$ norm of $P_H(A,B)$ and its Moore-Penrose Inverse $P^+_H(A,B)$ can be written as:

\begin{align}
&\|P_H(A,B)\|_2 = \left(\sigma_{max}(W^u_H(A,B))\right)^{\tfrac{1}{2}}  && \|P^+_H(A,B)\|_2 = \left(\sigma_{min}(W^u_H(A,B))\right)^{-\tfrac{1}{2}}
\end{align}

Because we assume a sufficient degree of controllability, $P_H(A,B)$ is full row-rank. This implies that

\begin{align}
&\|P_H(A,B)\|_2 = \sqrt{\lambda_{max}(P_H(A,B)P^\top_H(A,B))} = \sqrt{ \sigma_{max}(W^u_H(A,B))}\\
&\left(\|P^+_H(A,B)\|_2\right)^{-1}  = \sqrt{\lambda_{min}(P_H(A,B)P^\top_H(A,B))} = \sqrt{ \sigma_{min}(W^u_H(A,B))}
\end{align}

Let $(A,B)$ be a fixed pair of $H$-controllable system matrices, and let $\F(A,B)$ denote the matrix

\begin{align}
\F(A,B) = \begin{bmatrix}\bm{C}& 0\\ 0 & \bm{D}\end{bmatrix}\begin{bmatrix}\bm{G}_u(A,B) \\ I \end{bmatrix}(I-P^+_H(A,B)P_H(A,B)).
\end{align}

Then, $\|\F^+(A,B)\|_2 \leq \sigma^{-1}_{min}(D)$.

For an arbitrary matrix $M$, $(\|M^+\|_2)^{-1}$ is equal to the smallest non-zero singular eigenvalue of $M$ (we will denote this quantity as $\sigma_{-1}(M)$). Thus, in order to bound $\|M^+\|_2$ from above, we have to bound $\sigma_{-1}(M)$ from below. Denote $\bm{L}$ as the matrix
$$ \bm{L}:=\begin{bmatrix}\bm{C}& 0\\ 0 & \bm{D}\end{bmatrix}\begin{bmatrix}\bm{G}_u(A,B) \\ I \end{bmatrix} $$
and notice that it is full column rank and has rank of $H\times n_u$. The projection $\Pi_{\cl{N}(P_H)}:=(I-P^+_H(A,B)P_H(A,B))$ has rank $H\times n_u-n_x$ due the assumption of $H$-controllability. Hence, $\bm{F}=\bm{L}\Pi_{\cl{N}(P_H)}$ is full column rank with rank $r_{\F} := H\times n_u-n_x$ and has a null space $\cl{N}(\F)$ of dimension $n_x$. From these observations, we can equivalently say that $\sigma_{-1}(\F)$ is the $r_{\F}$th largest (or equivalently $n_x+1$ smallest) singular eigenvalue of $\F$. Using the Minimax principle, we can therefore write:

\begin{align}
\sigma_{-1}(\F) &= \max\limits_{\mathrm{proj. }\Pi, \text{ s.t.: } \mathrm{rank}(\Pi) = r_{\F}} \min\limits_{x\text{ s.t.: }\|\Pi x\|=1} x^\top \Pi \F^\top \F \Pi x \\
&= \max\limits_{\mathrm{proj. }\Pi, \text{ s.t.: } \mathrm{rank}(\Pi) = r_{\F}} \min\limits_{x\text{ s.t.: }\|\Pi x\|=1} x^\top \Pi \Pi_{\cl{N}(P_H)} \bm{L}^\top \bm{L} \Pi_{\cl{N}(P_H)}  \Pi x
\end{align}

Now recall that $\Pi_{\cl{N}(P_H)}$ is of rank $r_{\F}$, hence it is a feasible choice for the variable $\Pi$ of the outer optimization problem. This leads to the bound

\begin{align}
\sigma_{-1}(\F) &\geq \min\limits_{x\text{ s.t.: }\|\Pi_{\cl{N}(P_H)} x\|=1} x^\top \Pi_{\cl{N}(P_H)} \bm{L}^\top \bm{L} \Pi_{\cl{N}(P_H)} x \\
&\geq \min\limits_{z\text{ s.t.: }\|z\|=1} z^\top \bm{L}^\top \bm{L} z = \sigma_{min}(\bm{L})
\end{align}

We obtain a simple, but possibly conservative, lower bound on $\sigma_{min}(\bm{L})$ as follows:

\begin{align*}
&\sigma^2_{min}(\bm{L}) = \min\limits_{z\text{ s.t.: }\|z\|=1} \|\bm{L} z\|^2_2 = \min\limits_{z\text{ s.t.: }\|z\|=1} \|\bm{C}\bm{G}_u(A,B) z\|^2_2 + \|\bm{D}z\|^2_2 \geq \sigma^2_{min}(\bm{C}\bm{G}_u(A,B)) + \sigma^2_{min}(\bm{D})\\
&\implies \quad \sigma_{min}(\bm{L}) \geq \sigma_{min}(\bm{D})
\end{align*}

Finally, this provides us with the final result: $\|\F^+(A,B)\|_2 = \sigma^{-1}_{-1}(\bm{F}) \leq  \sigma^{-1}_{min}(\bm{L}) \leq \sigma^{-1}_{min}(\bm{D})$

We obtain an upper bound for $\|\g\|_2$, as a corollary of the previous three Lemmas:

Let $(A,B)$ be a fixed pair of $H$-controllable system matrices. Let $\g = \bm{L}\bm{u}^*_c$, where $\bm{L}$ and $\bm{u}^*_{c}$ are defined as:

\begin{align}
& \bm{L}:=\begin{bmatrix}\bm{C}& 0\\ 0 & \bm{D}\end{bmatrix}\begin{bmatrix}\bm{G}_u(A,B) \\ I \end{bmatrix} & &\bm{u}^*_{c} :=  P^{+}_H A^H e_j = P^\top_HW^{-1}_HA^H e_j.
\end{align}

Then, it holds:

$$\|\bm{g}\|_2 \leq \left(\|C\|_2 \sqrt{H}\sigma^{\tfrac{1}{2}}_{max}(W^u_H) + \|D\|_2\right)\sigma^{-\tfrac{1}{2}}_{min}(W^u_H)\alpha_H $$

where $\alpha_H := \max_{0\leq k\leq H} \|A^k\|_2$

§.§ The final bound

With the results of the last section, we can now bound the constants $\Gamma_A$ and $\Gamma_B$ used in thm:mainsensitivity. Rather than writing the explicit form of the constants we shall only analyze how they scale with system parameters. $\Gamma_A$, $\Gamma_B$ are defined as

\begin{align*}
\Gamma_A &= \kappa_{CD}\Gamma'_1 + \kappa_{CD}\Gamma'_2 \|B_1\|_2\|\bm{G}_w(A_1)\|_2 \|\bm{G}_w(A_2)\|_2 \\
\Gamma_B &= \kappa_{CD}\Gamma'_2 \|\bm{G}_w(A_2)\|_2,
\end{align*}

where $\Gamma'_1$, $\Gamma'_2$ are dominated by the terms:

\begin{align*}
\Gamma'_1 &\sim \cl{O}\left(\alpha_{H,1} \alpha_{H,2}H\left\|\bm{G}_{u,2}\right\|_2\|P^+_{H,2}\|_2\right)\\
\Gamma'_2 &\sim \cl{O}\left(\|\g\|_2(\|\F^+_1\|_2 + \|\F^+_2\|_2)\|P^+_{H,1}\|_2 \|P^+_{H,2}\|_2(\|P_{H,1}\|_2 + \|P_{H,2}\|_2)(1+\|\bm{G}_{u,1}\|_2)\right)
\end{align*}

Recall the FIR horizon $H$ and the constants $ \underline{\sigma}_u,  \overline{\sigma}_u$,  $\underline{\sigma}_w$,  $\overline{\sigma}_w$, $\kappa_{CD}$ and revisit the collection of bounds we have derived:

* $\|\bm{G}_{u}(A,B)\|_2 \leq \sqrt{H \sigma_{max}(W^u_H(A,B))}$,  $\|\bm{G}_{w}(A)\|_2 \leq \sqrt{H \sigma_{max}(W^w_H(A))} $

* $\|P_H(A,B)\|_2 = \left(\sigma_{max}(W^u_H(A,B))\right)^{\tfrac{1}{2}},\quad \|P^+_H(A,B)\|_2 = \left(\sigma_{min}(W^u_H(A,B))\right)^{-\tfrac{1}{2}}$

* $\|\F^+(A,B)\|_2 \leq \sigma^{-1}_{min}(D)$

* $\|\bm{g}\|_2 \leq \left(\|C\|_2 \sqrt{H}\sigma^{\tfrac{1}{2}}_{max}(W^u_H) + \|D\|_2\right)\sigma^{-\tfrac{1}{2}}_{min}(W^u_H)\alpha_H $

* $\alpha_H := \max_{0\leq k\leq H} \|A^k\|_2$

* $\kappa_{CD} = \frac{\max\{\sigma_{max}(C),\sigma_{max}(D)\}}{\min\{\sigma_{min}(C),\sigma_{min}(D)\}}$

Per lem:baseCtrl, we can bound the singular values of the arbitrary grammians using $ \underline{\sigma}_u,  \overline{\sigma}_u$,  $\underline{\sigma}_w$,  $\overline{\sigma}_w$:

\begin{align}
\underline{\sigma}_u \leq &\sigma_{min}(W^u_H(A_1,B_1)), \sigma_{min}(W^u_H(A_2,B_2)),\\
&\sigma_{max}(W^u_H(A_1,B_1)), \sigma_{max}(W^u_H(A_2,B_2)) \leq \overline{\sigma}_u \\
\underline{\sigma}_w \leq &\sigma_{min}(W^w_H(A_1)), \sigma_{min}(W^w_H(A_2)),\\
&\sigma_{max}(W^w_H(A_1)), \sigma_{max}(W^w_H(A_2)) \leq \overline{\sigma}_w
\end{align}

Then, we obtain

\begin{align}
&\Gamma'_2 = \cl{O}\left(
\alpha_H \:\kappa_{CD}\: H \:\left(\frac{\overline{\sigma}_u}{\underline{\sigma}_u}\right)^{\frac{3}{2}}\right) && \Gamma'_1 = \cl{O}\left( \alpha^2_H H^{\frac{3}{2}}\left(\frac{\overline{\sigma}_u}{\underline{\sigma}_u}\right)^{\frac{1}{2}} \right)
\end{align}

and finally

\begin{align}
&\Gamma_A = \cl{O}\left(
\alpha^2_H \:\kappa^2_{CD}\:\|B_1\|_2\: H^2 \:\left(\frac{\overline{\sigma}_u}{\underline{\sigma}_u}\right)^{\frac{3}{2}} \overline{\sigma}_w\right) && \Gamma_B = \cl{O}\left(
\alpha_H \:\kappa^2_{CD}\: H^{\frac{3}{2}} \:\left(\frac{\overline{\sigma}_u}{\underline{\sigma}_u}\right)^{\frac{3}{2}} \overline{\sigma}_w^{\frac{1}{2}}\right)
\end{align}

Let $C,D \succ 0$, and let $\cl{S}$ be a compact set of controllable systems with known FIR horizon $H$ and constants $ \underline{\sigma}_u,  \overline{\sigma}_u$,  $\underline{\sigma}_w$, $\overline{\sigma}_w$ as defined in lem:baseCtrl. Then there are fixed constants $\Gamma_A,\Gamma_B$, such that for any two pairs of system matrices $(A_1,B_1), (A_2,B_2) \in \cl{S}$ the corresponding $\cl{H}_2$ optimal SLS-solutions of problem $S_j$ ($j$ arbitrary), denoted $\phi^{*j}_1$ and $\phi^{*j}_2$, satisfy the following inquality:

\begin{align}
\|\phi^{j*}_1-\phi^{j*}_2\|_2 \leq \Gamma_A\|A_1-A_2\|_F + \Gamma_B \|B_1-B_2\|_F.
\end{align}

Furthermore, $\Gamma_A$ and $\Gamma_B$ satisfy

\begin{align}
&\Gamma_A = \cl{O}\left(
\alpha^2_H \:\kappa^2_{CD}\:\beta\: H^2 \:\left(\frac{\overline{\sigma}_u}{\underline{\sigma}_u}\right)^{\frac{3}{2}} \overline{\sigma}_w\right) && \Gamma_B = \cl{O}\left(
\alpha_H \:\kappa^2_{CD}\: H^{\frac{3}{2}} \:\left(\frac{\overline{\sigma}_u}{\underline{\sigma}_u}\right)^{\frac{3}{2}} \overline{\sigma}_w^{\frac{1}{2}}\right),
\end{align}

where $\beta:= \max\limits_{(A,B)\in\cl{S}}\|B\|_2$ and $\kappa_{CD}$ stands for

$$\kappa_{CD} = \frac{\max\{\sigma_{max}(C),\sigma_{max}(D)\}}{\min\{\sigma_{min}(C),\sigma_{min}(D)\}}$$

§ VECTOR SUBSYSTEMS

For ease of exposition, Algorithm <ref> was stated for the case when each subsystem $i \in [N]$ in the network has scalar dynamics. Below we describe how Algorithm <ref> can be implemented for general vector subsystems that are heterogeneous.
Recall that we consider a network of $N$ subsystems where each subsystem $i$ has $x^i(t) \in \mathbb{R}^{n_i}$ and $u^i(t) \in \mathbb{R}^{m_i}$.

Communication and locality constraints.
The communication graph and the communication matrix $\mathcal{C}\in \{0,1\}^{N \times N}$ can be identically defined for these vector subsystems as the scalar case in section <ref>, e.g., having vector subsystems does not change the communication graph and all definitions according to the communication graph.  Similarly, the vector subsystems do not affect our definition of the locality constraints as described in <Ref>.

When we write $\Phi^x[k] \in \mathcal{C}$, where $\Phi^x[k] \in \R^{n_x \times n_x}$ and $\mathcal{C}\in \{0,1\}^{N\times N}$ with $n_x = \sum_{i=1}^N n_i$, we mean that the $i,j^{\text{th}}$ block matrix of $\Phi^x[k]$ (of size $n_i$ by $n_j$) is nonzero if and only if $\mathcal{C}(i,j) =1$.

If the system is not fully actuated where $B\in \R^{n_x \times n_u}$ with $n_u < n_x$, then $\Phi^u[k] \in \R^{n_u \times n_x}$ and thus is a wide matrix.
In this case, we write $\Phi^u[k] \in \mathcal{C}$ to mean that each $i^{\text{th}}$ row of $\Phi^u[k] $ (corresponding to the actuated position $u^i(t)$) conforms to the sparsity pattern of the $i^{\text{th}}$ row of $\mathcal{C}$, where if $u^i(t)$ belongs to the $\ell^{\text{th}}$ subsystem, $\Phi^u[k](i,k) \not=  0$ if and only if $\mathcal{C}(\ell,q) =1$ with $x^k(t)$ belonging to subsystem $q$.

The CONSIST component (Algorithm <ref>) of Algorithm <ref> does not change for the vector case, because (<ref>) is general for any local state dimension. Instead of computing each row of $A,B$, each subsystem now computes multiple rows of $A,B$ that corresponds to its local dynamics in (<ref>).

The CONTROL component (Algorithm <ref>) of Algorithm <ref> is where modification is required.
Before we state how to generalize CONTROL to the vector case, we need a few more definitions.
Specifically, define index set $\{i_k\}$ for $i\in[N]$ and $k\in[n_i]$ as the state index for subsystem $i$, where each $i_k$ is the position of the $k^\text{th}$ component of $x^i$ in the global state vector $x$. A control index can be analogously defined for subsystem $i$.

In the first step of CONTROL where each subsystem synthesizes a single closed-loop column in the scalar case, each subsystem now synthesizes multiple closed-loop columns. In particular, subsystem $i$ computes all $i_k^{\text{th}}$ closed-loop columns using (<ref>). During communication, all $i_k$ columns from subsystem $i$ is transmitted to its outgoing communication neighbors.

In the second step of CONTROL where each subsystem computes its local control action, subsystem $i$ now computes $\hat{w}^{i_k}(t)$ and $u^{i_k}(t)$ for all its state and control index in (<ref>), rather than the single index $i$ in the scalar case.
This also means that subsystem $i$ uses multiple closed-loop columns transmitted from other subsystems.
In terms of the distributed SLS controller implementation (described in <Ref>), this means that instead of performing a single row-vector multiplication during (<ref>), subsystem $i$ now performs matrix-vector multiplication (equivalently, multiple row-vector multiplication) over each of its state and control index $i_k$. In particular, (<ref>) now requires the closed-loop columns from subsystem $j$ over all state and control index $\{j_k\}$ of subsystem $j$ for $j \in \din{i}$.

Stability guarantee.
It is straightforward to verify that the analysis in <Ref> can be carried through in near identical fashion. The effect of general vector subsystem  manifests in the change of the constants such as $\bar{n}$, where the constants are larger than the scalar case. In general the size of each subsystem is still much smaller than the dimension of the global system.