# Simultaneous Input and State Set-Valued $\mathcal{H}_{\infty}$-Observers For
Linear Parameter-Varying Systems

Mohammad Khajenejad Sze Zheng Yong
M. Khajenejad and S.Z. Yong are with the School for Engineering of Matter,
Transport and Energy, Arizona State University, Tempe, AZ, USA (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>).

###### Abstract

A fixed-order set-valued observer is presented for linear parameter-varying
systems with bounded-norm noise and under completely unknown attack signals,
which simultaneously finds bounded sets of states and unknown inputs that
include the true state and inputs. The proposed observer can be designed using
semidefinite programming with LMI constraints and is optimal in the minimum
$\mathcal{H}_{\infty}$-norm sense. We show that the strong detectability of
each _constituent_ linear time-invariant system is a necessary condition for
the existence of such an observer, as well as the boundedness of set-valued
estimates. Furthermore, sufficient conditions are provided for the upper
bounds of the estimation errors to converge to steady state values and
finally, the results of such a set-valued observer are exhibited through an
illustrative example.

##  I Introduction

The security of Cyber-Physical Systems (CPS) is emerging as an extremely
critical and important issue. Since physical and software components are
deeply intertwined in CPS, such systems are potentially vulnerable to
adversarial attacks, which could be harmful for both the physical systems and
their operators. Given that adversarial attackers may behave strategically,
there are many potential avenues through which they can cause harm, steal
information/power, etc.

Misleading the system operator by inserting counterfeit data into sensor and
actuator signals (false data injection) is among the most common and extensive
attacks on CPS. Hence, a significant amount of effort has been invested in new
designs for estimation and control against false data injection attacks. Due
to the nature of the attack signals, it is not justifiable to impose any kind
of restrictive assumptions on them (e.g., stochastic with normal distribution
or deterministic with bounded norm) i.e., they can be _anything_ , so none of
the classical Kalman filtering based methods are applicable. Moreover,
considering the complicated structure of CPS (several physical, computing and
communication components), they can be modeled and represented more
realistically by time-varying and nonlinear dynamic systems rather than linear
time-invariant ones. Taking all these facts into account, this paper attempts
to design a resilient observer for a particular class of linear time-varying
systems known as linear parameter-varying systems, which to the best of our
knowledge is a novel approach.
_Literature review._ There are several different frameworks for simultaneous
input and state estimation of linear time-varying stochastic systems with
unknown inputs, assuming that the noise signals are Gaussian and white. The
authors in [1, 2, 3, 4] apply Kalman filtering inspired recursive filter
design approaches (modified versions of unbiased minimum-variance estimation
methods), where some additional assumptions are needed to guarantee the
stability of the filters, while [5] uses a modified double-model adaptive
estimation method. However, these Kalman filtering inspired approaches are not
applicable for set-membership estimation problems (cf. [6] for a comprehensive
discussion), as is considered in this paper.

In the context of attack-resilient estimation, where adversarial signals can
be malicious and strategic and thus, bounds on the unknown inputs/attack
signals are completely unknown, there have been a number of proposed
approaches in the literature for systems with bounded errors (e.g., [7, 8,
9]), but all of them only consider point estimates, i.e, the most likely or
best single estimate as opposed to a set-valued estimate. Specifically, the
work in [7] only computes error bounds for the initial state and [8] assumes
zero initial states and does not consider any optimality criteria. The author
in [6] and references therein discussed the advantages of set-valued observers
(when compared to point estimators) in terms of providing hard accuracy
bounds, which are important to guarantee safety [10]. In addition, the use of
fixed-order set-valued methods can help decrease the complexity of optimal
observers[11], which grows with time. Hence, the work in [6] presents a fixed-
order set-valued observer for linear time-invariant discrete time systems with
bounded errors, that simultaneously finds bounded sets of compatible states
and unknown inputs that are optimal in the minimum $\mathcal{H}_{\infty}$-norm
sense, i.e., with minimum average power amplification, which we aim to
generalize in this paper for linear parameter-varying systems.
_Contributions._ We propose a novel fixed-order set-valued observer for linear
parameter-varying systems with unknown input and bounded noise signals that
simultaneously finds bounded sets of states and unknown inputs that contain
the true state and unknown input and are compatible/consistent with the
measurement outputs. Specifically, we consider linear parameter-varying system
dynamics that can be presented as a convex combination of linear time-
invariant _constituent_ dynamics. In addition, we provide necessary conditions
for the boundedness of the set-valued estimates. We further prove the
optimality of the filter in the minimum $\mathcal{H}_{\infty}$-norm sense,
i.e., minimum average power amplification, by converting the corresponding
problem into a tractable formulation using semi-definite programming with LMI
constraints that is readily implementable using off-the-shelf optimization
solvers. We also show that strong detectability of each constituent system is
a necessary condition for the existence of such an
$\mathcal{H}_{\infty}$-observer. Then, we provide some sufficient conditions
for the convergence of upper bounds of the state and input estimation errors
to steady state and for obtaining these steady state bounds. Finally, we
demonstrate the effectiveness of our proposed set-valued observer through an
illustrative example.

_Notation._ $\mathbb{R}^{n}$ denotes the $n$-dimensional Euclidean space and
$\mathbb{N}$ nonnegative integers. For a vector $v\in\mathbb{R}^{n}$ and a
matrix $M\in\mathbb{R}^{p\times q}$, $\|v\|\triangleq\sqrt{v^{\top}v}$ and
$\|M\|$ denote their (induced) 2-norm. Moreover, the transpose, inverse,
Moore-Penrose pseudoinverse and rank of $M$ are given by $M^{\top}$, $M^{-1}$,
$M^{\dagger}$ and ${\rm rk}(M)$. For a symmetric matrix $S$, $S\succ 0$
($S\succeq 0$) is positive (semi-)definite.

##  II Problem Statement

_System Assumptions._ Consider the following linear parameter-varying
discrete-time bounded-error system:

$\displaystyle\begin{split}x_{k+1}&=\textstyle\sum_{i=1}^{N}\lambda_{i,k}({A}^{i}x_{k}+{B}^{i}u_{k}+{w}^{i}_{k})+Gd_{k},\\\
y_{k}&={C}x_{k}+\textstyle\sum_{i=1}^{N}\lambda_{i,k}({D}^{i}u_{k}+{v}^{i}_{k})+Hd_{k},\end{split}$
(1)

where $\lambda_{i,k}$ is known and satisfies $0\leq\lambda_{i,k}\leq
1,\sum_{i=1}^{N}\lambda_{i,k}=1,\forall k$. $x_{k}\in\mathbb{R}^{n}$ is the
state vector at time $k\in\mathbb{N}$, $u_{k}\in\mathbb{R}^{m}$ is a known
input vector, $d_{k}\in\mathbb{R}^{p}$ is an unknown input vector, and
$y_{k}\in\mathbb{R}^{l}$ is the measurement vector. The process noise
${w}^{i}_{k}\in\mathbb{R}^{n}$ and the measurement noise
${v}^{i}_{k}\in\mathbb{R}^{l}$ are assumed to be bounded and $\ell_{\infty}$
sequences, with $\|{w}^{i}_{k}\|\leq\eta_{w}$ and
$\|{v}^{i}_{k}\|\leq\eta_{v}$.We also assume an estimate $\hat{x}_{0}$ of the
initial state $x_{0}$ is available, where
$\|\hat{x}_{0}-x_{0}\|\leq\delta_{0}^{x}$. The matrices ${A}^{i}$, ${B}^{i}$,
${C}$, ${D}^{i}$, $G$ and $H$ are known for $i\in\left\\{1,2,\dots,N\right\\}$
and of appropriate dimensions, where $G$ and $H$ are matrices that encode the
_locations_ through which the unknown input or attack signal can affect the
system dynamics and measurements and $N$ is the number of _constituent_
systems. Note that no assumption is made on $H$ to be either the zero matrix
(no direct feedthrough), or to have full column rank when there is direct
feedthrough. Without loss of generality, we assume that ${\rm
rk}[G^{\top}\;H^{\top}]=p$, $n\geq l\geq 1$, $l\geq p\geq 0$, $m\geq 0$ and
each $({A}^{i},{B}^{i},{C},{D}^{i},G,H),i\in\left\\{1,2,\dots,N\right\\}$
represents a linear time-invariant constituent system:

$\displaystyle\begin{array}[]{ll}x^{i}_{k+1}&={A}^{i}x^{i}_{k}+{B}^{i}u_{k}+Gd_{k}+{w}^{i}_{k},\\\
y^{i}_{k}&={C}x_{k}+{D}^{i}u_{k}+Hd_{k}+{v}^{i}_{k}.\end{array}$ (4)

_Unknown Input (or Attack) Signal Assumptions._ The unknown inputs $d_{k}$ are
not constrained to be a signal of any type (random or strategic) nor to follow
any model, thus no prior ‘useful’ knowledge of the dynamics of $d_{k}$ is
available (independent of $\\{d_{\ell}\\}$ $\forall k\neq\ell$,
$\\{w_{\ell}\\}$ and $\\{v_{\ell}\\}$ $\forall\ell$). We also do not assume
that $d_{k}$ is bounded or has known bounds and thus, $d_{k}$ is suitable for
representing adversarial attack signals.

The simultaneous input and state set-valued observer design problem can be
stated as follows:

######  Problem 1.

Given a linear parameter-varying discrete-time bounded-error system with
unknown inputs (1) , design an optimal and stable filter that simultaneously
finds bounded sets of compatible states and unknown inputs in the minimum
$\mathcal{H}_{\infty}$-norm sense, i.e., with minimum average power
amplification.

##  III Preliminary Material

###  III-A System Transformation

In order to decouple the output equation into two components, first a
transformation is carried out for each of the constituent subsystems, one with
a full rank direct feedthrough matrix and the other without direct
feedthrough. Note that this similarity transformation is similar to the one in
[6] and is not the same as the one in [4], which is no longer applicable as it
was based on the noise error covariance.

Let $p_{H}\triangleq{\rm rk}(H)$. Using singular value decomposition, we
rewrite the direct feedthrough matrix $H$ as
$H=\begin{bmatrix}U_{1}&U_{2}\end{bmatrix}\begin{bmatrix}\Sigma&0\\\
0&0\end{bmatrix}\begin{bmatrix}V_{1}^{\,\top}\\\ V_{2}^{\,\top}\end{bmatrix}$,
where $\Sigma\in\mathbb{R}^{p_{H}\times p_{H}}$ is a diagonal matrix of full
rank, $U_{1}\in\mathbb{R}^{l\times p_{H}}$,
$U_{2}\in\mathbb{R}^{l\times(l-p_{H})}$, $V_{1}\in\mathbb{R}^{p\times p_{H}}$
and $V_{2}\in\mathbb{R}^{p\times(p-p_{H})}$, while
$U\triangleq\begin{bmatrix}U_{1}&U_{2}\end{bmatrix}$ and
$V\triangleq\begin{bmatrix}V_{1}&V_{2}\end{bmatrix}$ are unitary matrices.
When there is no direct feedthrough, $\Sigma$, $U_{1}$ and $V_{1}$ are empty
matricesaaa Based on the convention that the inverse of an empty matrix is an
empty matrix and the assumption that operations with empty matrices are
possible., and $U_{2}$ and $V_{2}$ are arbitrary unitary matrices.

Then, we decouple the unknown input into two orthogonal components:

$\displaystyle d_{1,k}=V_{1}^{\top}d_{k},\quad d_{2,k}=V_{2}^{\top}d_{k}.$ (5)

Considering that $V$ is unitary, $d_{k}=V_{1}d_{1,k}+V_{2}d_{2,k}$ and we can
represent the system (1) as:

$\displaystyle x_{k+1}$
$\displaystyle=\textstyle\sum_{i=1}^{N}\lambda_{i,k}({A}^{i}x_{k}+{B}^{i}u_{k}+{w}^{i}_{k})+G_{1}d_{1,k}+G_{2}d_{2,k},$
$\displaystyle y_{k}$
$\displaystyle={C}x_{k}+\textstyle\sum_{i=1}^{N}\lambda_{i,k}({D}^{i}u_{k}+{v}^{i}_{k})+H_{1}d_{1,k}$
(6)

where $G_{1}\triangleq GV_{1}$, $G_{2}\triangleq GV_{2}$ and $H_{1}\triangleq
HV_{1}=U_{1}\Sigma$. Next, the output $y_{k}$ is decoupled using a nonsingular
transformation
$T=\begin{bmatrix}T_{1}^{\top}&T_{2}^{\top}\end{bmatrix}^{\top}\triangleq
U^{\top}=\begin{bmatrix}U_{1}&U_{2}\end{bmatrix}^{\top}$ to get
$z_{1,k}\in\mathbb{R}^{p_{H}}$ and $z_{2,k}\in\mathbb{R}^{l-p_{H}}$ given by

$\displaystyle\begin{array}[]{lll}z_{1,k}&\triangleq
T_{1}y_{k}=U_{1}^{\top}y_{k}\\\ &={C}_{1}x_{k}+\Sigma
d_{1,k}+\sum_{i=1}^{N}\lambda_{i,k}{D}^{i}_{1}u_{k}+\sum_{i=1}^{N}\lambda_{i,k}{v}^{i}_{1,k}\\\
z_{2,k}&\triangleq T_{2}y_{k}=U_{2}^{\top}y_{k}\\\
&={C}_{2}x_{k}+\sum_{i=1}^{N}\lambda_{i,k}{D}^{i}_{2}u_{k}+\sum_{i=1}^{N}\lambda_{i,k}{v}^{i}_{2,k}\vspace{-0.5cm}\end{array}$
(11)

where ${C}_{1}\triangleq U_{1}^{\top}{C}$, ${C}_{2}\triangleq
U_{2}^{\top}{C}$, ${D}^{i}_{1}\triangleq U_{1}^{\top}{D}^{i}$,
${D}^{i}_{2}\triangleq U_{2}^{\top}{D}^{i}$, ${v}^{i}_{1,k}\triangleq
U_{1}^{\top}{v}^{i}_{k}$ and ${v}^{i}_{2,k}\triangleq
U_{2}^{\top}{v}^{i}_{k}$. This transform is also chosen such that
$\|\begin{bmatrix}{{v}^{i}_{1,k}}^{\top}&{{v}^{i}_{2,k}}^{\top}\end{bmatrix}^{\top}\|=\|U^{\top}{v}^{i}_{k}\|=\|{v}^{i}_{k}\|$.

##  IV Fixed-Order Simultaneous Input and State Set-Valued Observers

###  IV-A Set-Valued Observer Design

We consider a recursive three-step set-valued observer design. This design
utilizes a similar framework as in [6] and contains an _unknown input
estimation_ step that uses the current measurement and the set of compatible
states to estimate the set of compatible unknown inputs, a _time update_ step
which propagates the compatible set of states based on the system dynamics,
and a _measurement update_ step that uses the current measurement to update
the set of compatible states. To sum up, our target is to design a three-step
recursive set-valued observer of the form:

Unknown Input Estimation: $\displaystyle\
\hat{D}_{k-1}=\mathcal{F}_{d}(\hat{X}_{k-1},u_{k}),$ Time Update:
$\displaystyle\quad\
\hat{X}^{\star}_{k}=\mathcal{F}_{x}^{\star}(\hat{X}_{k-1},\hat{D}_{k-1},u_{k}),$
Measurement Update: $\displaystyle\quad\
\hat{X}_{k}=\mathcal{F}_{x}(\hat{X}^{\star}_{k},u_{k},y_{k}),$

where $\mathcal{F}_{d}$, $\mathcal{F}^{\star}_{x}$ and $\mathcal{F}_{x}$ are
to-be-designed set mappings, while $\hat{D}_{k-1}$, $\hat{X}^{\star}_{k}$ and
$\hat{X}_{k}$ are the sets of compatible unknown inputs at time $k-1$,
propagated, and updated states at time $k$, correspondingly. It is important
to note that $d_{2,k}$ cannot be estimated from $y_{k}$ since it does not
affect $z_{1,k}$ and $z_{2,k}$. Thus, the only estimate we can obtain in light
of (11) is a (one-step) delayed estimate of $\hat{D}_{k-1}$. The reader may
refer to a previous work [3] for a complete discussion on when a delay is
absent or when we can expect further delays. Similar to [12, 10, 6], as the
complexity of optimal observers increases with time, only the fixed-order
recursive filters will be considered. In particular, we choose set-valued
estimates of the form:

$\displaystyle\begin{array}[]{rl}\hat{D}_{k-1}&=\\{d\in\mathbb{R}^{p}:\|d_{k-1}-\hat{d}_{k-1}\|\leq\delta^{d}_{k-1}\\},\\\
\hat{X}^{\star}_{k}&=\\{x\in\mathbb{R}^{n}:\|x_{k}-\hat{x}^{\star}_{k|k}\|\leq\delta^{x,\star}_{k}\\},\\\
\hat{X}_{k}&=\\{x\in\mathbb{R}^{n}:\|x_{k}-\hat{x}_{k|k}\|\leq\delta^{x}_{k}\\}.\end{array}$

In other words, we restrict the estimation errors to balls of norm $\delta$.
In this setting, the observer design problem is equivalent to finding the
centroids $\hat{d}_{k-1}$, $\hat{x}^{\star}_{k|k}$ and $\hat{x}_{k|k}$ as well
as the radii $\delta^{d}_{k-1}$, $\delta_{k}^{x,\star}$ and $\delta_{k}^{x}$
of the sets $\hat{D}_{k-1}$, $\hat{X}^{\star}_{k}$ and $\hat{X}_{k}$,
respectively. In addition, we limit our attention to observers for the
centroids $\hat{d}_{k-1}$, $\hat{x}^{\star}_{k|k}$ and $\hat{x}_{k|k}$ that
belong to the class of three-step recursive filters given in [2] and [4],
defined as follows for each time $k$ (with $\hat{x}_{0|0}=\hat{x}_{0}$):

_Unknown Input Estimation_ :

$\displaystyle\hat{d}_{1,k}=M_{1}(z_{1,k}-{C}_{1}\hat{x}_{k|k}-\textstyle\sum_{i=1}^{N}\lambda_{i,k}{D}^{i}_{1}u_{k}),$
(12)
$\displaystyle\hat{d}_{2,k-1}=M_{2}(z_{2,k}-{C}_{2}\hat{x}_{k|k-1}-\textstyle\sum_{i=1}^{N}\lambda_{i,k}{D}^{i}_{2}u_{k}),$
(13) $\displaystyle\hat{d}_{k-1}=V_{1}\hat{d}_{1,k-1}+V_{2}\hat{d}_{2,k-1}.$
(14)

_Time Update_ :

$\displaystyle\hat{x}_{k|k-1}$
$\displaystyle=\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}({A}^{i}\hat{x}_{k-1|k-1}+{B}^{i}u_{k-1})+G_{1}\hat{d}_{1,k-1},$
(15) $\displaystyle\hat{x}^{\star}_{k|k}$
$\displaystyle=\hat{x}_{k|k-1}+G_{2}\hat{d}_{2,k-1}.$ (16)

_Measurement Update_ :

$\displaystyle\begin{array}[]{rl}\hat{x}_{k|k}&=\hat{x}^{\star}_{k|k}+L(y_{k}-{C}\hat{x}^{\star}_{k|k}-\textstyle\sum_{i=1}^{N}\lambda_{i,k}{D}^{i}u_{k})\\\
&=\hat{x}^{\star}_{k|k}+\tilde{L}(z_{2,k}-{C}_{2}\hat{x}^{\star}_{k|k}-\textstyle\sum_{i=1}^{N}\lambda_{i,k}{D}^{i}_{2}u_{k}),\end{array}$
(19)

where $L\in\mathbb{R}^{n\times l}$, $\tilde{L}\triangleq
LU_{2}\in\mathbb{R}^{n\times(l-p_{H})}$, $M_{1}\in\mathbb{R}^{p_{H}\times
p_{H}}$ and $M_{2}\in\mathbb{R}^{(p-p_{H})\times(l-p_{H})}$ are observer gain
matrices that are designed according to Theorem 1. The main result in Theorem
1 is derived by minimizing the “volume” of the set of compatible states and
unknown inputs, quantified by the radii $\delta^{d}_{k-1}$,
$\delta_{k}^{x,\star}$ and $\delta_{k}^{x}$. Note also that we applied
$L=LU_{2}U_{2}^{\top}=\tilde{L}U_{2}^{\top}$ from Lemma 1 into (19). The state
and input estimation errors are defined as $\tilde{x}_{k|k}\triangleq
x_{k}-\hat{x}_{k|k},\tilde{d}_{k-1}\triangleq
d_{k-1}-\hat{d}_{k-1},\tilde{d}_{1,k-1}\triangleq
d_{1,k-1}-\hat{d}_{1,k-1},\tilde{d}_{2,k-1}\triangleq
d_{2,k-1}-\hat{d}_{2,k-1}$ respectively. In Lemmas 1 and 2, we will provide
necessary conditions for boundedness of estimation errors and sufficient
conditions for stability of the observer. All the proofs are provided in the
Appendix.

######  Lemma 1 (Necessary Conditions for Boundedness of Set-Valued Estimates
[6, Lemma 1]).

The input and state estimation errors, ($\tilde{d}_{k-1}$ and
$\tilde{x}_{k|k}$), are bounded for all $k$ (i.e., the set-valued estimates
are bounded with radii
$\delta^{d}_{k-1},\delta_{k}^{x,\star},\delta_{k}^{x}<\infty$), only if
$M_{1}\Sigma=I$, $p\leq l$, $M_{2}C_{2}G_{2}=I$ and $LU_{1}=0$ . Consequently,
${\rm rk}(C_{2}G_{2})=p-p_{H}$, $M_{1}=\Sigma^{-1}$,
$M_{2}=(C_{2}G_{2})^{\dagger}$ and
$L=LU_{2}U_{2}^{\top}=\tilde{L}U_{2}^{\top}$.

######  Lemma 2 (Sufficient Conditions for Observer Stability).

A sufficient condition for the stability of the set-valued observer is that
$(\overline{A}_{k},C_{2})$ is uniformly detectablebbb For conciseness, the
readers are referred to [13, Section 2] for the definition of uniform
detectability. A spectral test can be found in [14]. for each $k$, where
$\overline{A}_{k}\triangleq(I-G_{2}M_{2}C_{2})\hat{A}_{k}$ and
$\hat{A}_{k}\triangleq\sum_{i=1}^{N}\lambda_{i,k}{A}^{i}-G_{1}M_{1}C_{1}$.

###  IV-B Optimal $\mathcal{H}_{\infty}$-Observer

In this section, we provide sufficient conditions for the _existence_ of a
set-valued observer for system (1) with any sequence
$\\{\lambda_{i,k}\\}_{k=0}^{\infty}$ for all $i\in\\{1,2,\dots,N\\}$ that
satisfies $0\leq\lambda_{i,k}\leq 1,\sum_{i=1}^{N}\lambda_{i,k}=1,\forall k$
in the sense of $\mathcal{H}_{\infty}$ (i.e., minimizing the sum of squares of
the state estimation error sequence). Furthermore, we introduce a relatively
simple approach to find such an observer, which involves solving a semi-
definite program with Linear Matrix Inequalities (LMI) as constraints. We will
also show that given some structural conditions for the system, the upper
bounds of the estimation errors for both states and unknown inputs are
guaranteed to converge to steady state.

######  Theorem 1 ($\mathcal{H}_{\infty}$-Observer Design).

Suppose Lemma 1 holds and there exist matrices $Y$ and $S\succ 0$ with
appropriate dimensions such that

$\begin{bmatrix}S&(\overline{{A}^{i}})^{\top}(S-C_{2}^{\top}Y^{\top})&0&I\\\
*&S&\begin{bmatrix}S-YC_{2}&-Y\end{bmatrix}&0\\\ *&*&\eta I&0\\\ *&*&*&\eta
I\end{bmatrix}\succ 0$

for all $i\in\left\\{1,2,\dots,N\right\\}.$ Then, there exists an $\eta$
performance bounded $\mathcal{H}_{\infty}$-observer for system (1) with any
sequence $\\{\lambda_{i,k}\\}_{k=0}^{\infty}$ for all $i\in\\{1,2,\dots,N\\}$
that satisfies $0\leq\lambda_{i,k}\leq 1,\sum_{i=1}^{N}\lambda_{i,k}=1,\forall
k$ when using $\tilde{L}={S}^{-1}Y$, i.e.,
$\|T_{\tilde{x},w,v}\|\leq\eta^{2}$, where $T_{\tilde{x},w,v}$ is the transfer
function matrix that maps the noise signals
$\sum_{i=1}^{N}\lambda_{i,k}\begin{bmatrix}w_{k}^{i\,\top}&v_{k}^{i\,\top}\end{bmatrix}^{T}$
to the updated state estimation error $\tilde{x}_{k|k}\triangleq
x_{k}-\hat{x}_{k|k}$.

Furthermore, the optimal filter gain
$\tilde{L}={S^{\star}}^{-1}\tilde{Y}^{\star}$ with $\eta^{\star}$
$\mathcal{H}_{\infty}$-performance can be obtained from the following semi-
definite programming with LMI constraints:

$\displaystyle(\eta^{\star},S^{\star},Y^{\star})$
$\displaystyle\in\operatornamewithlimits{arg\ min}_{\eta,S,Y}\hskip
5.69046pt\eta$ $\displaystyle
s.t\begin{bmatrix}S&(\overline{{A}^{i}})^{\top}(S-C_{2}^{\top}Y^{\top})&0&I\\\
*&S&\begin{bmatrix}S-YC_{2}&-Y\end{bmatrix}&0\\\ *&*&\eta I&0\\\ *&*&*&\eta
I\end{bmatrix}\succ 0,$ $\displaystyle\hskip 106.69783pt\forall
i\in\left\\{1,2,..,N\right\\}.$ (20)

Although Theorem 1 equips us with an approach for designing an
$\mathcal{H}_{\infty}$-observer for the linear parameter-varying system in (1)
when one exists, it would still be valuable to find a _structural_ and
conveniently testable property for the constituent linear time-invariant
systems in (4) that is _necessary_ for the existence of such an observer.
Knowing such conditions would be beneficial in the sense that if they are
_not_ satisfied, the designer knows _a priori_ that there does not exist any
$\mathcal{H}_{\infty}$-observer for such an attacked system. This will be the
goal of Theorem 2.

######  Theorem 2 (Necessary Conditions for the Existence of an
$\mathcal{H}_{\infty}$-observer).

There exists a simultaneous state and unknown input
$\mathcal{H}_{\infty}$-observer for system (1) with any sequence
$\\{\lambda_{i,k}\\}_{k=0}^{\infty}$ for all $i\in\\{1,2,\dots,N\\}$ that
satisfies $0\leq\lambda_{i,k}\leq 1,\sum_{i=1}^{N}\lambda_{i,k}=1,\forall k$,
only if each $({A}^{i},G,{C},H)$ is strongly detectablecccFor brevity, the
readers may refer to [6] for the definition of strong detectability. for all
$i\in\left\\{1,2,\dots,N\right\\}$.

Next, we characterize the resulting radii $\delta^{x}_{k}$ and
$\delta^{d}_{k-1}$ when using the proposed $\mathcal{H}_{\infty}$-observer.

######  Theorem 3 (Radii of Set-Valued Estimates).

The radii $\delta^{x}_{k}$ and $\delta^{d}_{k-1}$ can be obtained as:

$\displaystyle\delta^{x}_{k}=\delta^{x}_{0}\theta^{k}+\overline{\eta}\textstyle\sum_{i=1}^{k}\theta^{i-1},$
$\displaystyle\delta^{d}_{k-1}=\beta\delta^{x}_{k-1}+\|V_{2}M_{2}C_{2}\|\eta_{w}+\big{[}\|(V_{2}M_{2}C_{2}G_{1}-V_{1})M_{1}T_{1}\|+\|V_{2}M_{2}T_{2}\|\big{]}\eta_{v},$

where
$\beta\triangleq\textstyle\max_{i\in\\{1,2,\dots,N\\}}\|V_{1}M_{1}C_{1}+V_{2}M_{2}C_{2}A_{e,i}\|,\Psi\triangleq
I-\tilde{L}C_{2},\Phi\triangleq
I-G_{2}M_{2}C_{2},A_{e,i}\triangleq\Psi\Phi(A^{i}-G_{1}M_{1}C_{1}),\\\
\theta\triangleq\max_{i\in\\{1,2,\dots,N\\}}\|A_{e,i}\|$.

The resulting fixed-order set-valued observer is summarized in Algorithm 1.

So far, we have designed an $\mathcal{H}_{\infty}$-observer for our linear
parameter-varying system and provided necessary conditions for the boundedness
of the set-valued estimates. It is worth mentioning that for the linear time-
invariant case in [6], strong detectability of the system is also a sufficient
condition for the convergence of the radii $\delta^{x}_{k}$ and
$\delta^{d}_{k-1}$ to steady state. In our parameter-varying case, even if all
constituent linear time-invariant systems are strongly detectable, there is no
guarantee that the radii converge. The reason is that the convergence hinges
on the stability of the product of _time-varying_ matrices (cf. proof of
Theorem 4), which is not guaranteed even if all the multiplicands are stable.
In the next theorem, we discuss some sufficient conditions for the convergence
of the radii to steady state.

######  Theorem 4 (Convergence).

Suppose the conditions of Theorem 1 hold. Then, the radii $\delta^{x}_{k}$ and
$\delta^{d}_{k-1}$ are convergent if $\|A_{e,i}\|<1$ for all
$i\in\left\\{1,2,\dots,N\right\\}$, where $A_{e,i}$ is defined in Theorem (3).
Moreover, the steady state radii is given by:

$\displaystyle\begin{array}[]{l}\displaystyle\lim_{k\to\infty}\delta^{x}_{k}=\frac{\overline{\eta}}{1-\theta},\\\
\displaystyle\lim_{k\to\infty}\delta^{d}_{k}=\frac{\overline{\eta}\beta}{1-\theta}+\eta_{w}\|V_{2}M_{2}C_{2}\|+\eta_{v}(\|V_{2}M_{2}T_{2}\|+\|R\|),\end{array}$

where
$\overline{\eta}\triangleq(\|\Gamma\|\eta_{v}+\|\Psi\Phi\|\eta_{w}),R\triangleq
V_{2}M_{2}C_{2}G_{1}M_{1}T_{1}-V_{1}M_{1}T_{1},\Gamma\triangleq-(\Psi\Phi
G_{1}M_{1}T_{1}+\Psi G_{2}M_{2}T_{2}+\tilde{L}T_{2})$.

######  Remark 1.

Alternatively, we can trade off between “optimality” of the observer and
“convergence” of the radii. We can iteratively find $\eta$ (e.g., by line
search) that satisfies the following feasibility problem:

$\displaystyle\begin{array}[]{rl}&{\rm Find}\hskip 5.69046pt(S,Y)\\\
&s.t\begin{bmatrix}S&*&0&I\\\
(S-YC_{2})\overline{A}^{i}&S&\begin{bmatrix}S-YC_{2}&-Y\end{bmatrix}&0\\\
*&*&\eta_{0}I&0\\\ *&*&*&\eta_{0}I\end{bmatrix}\succ 0,\forall
i\in\left\\{1,2,..,N\right\\},\end{array}$

as well as the sufficient condition in Theorem 4, i.e., $\|A_{e,i}\|<1$ for
all $i\in\left\\{1,2,\dots,N\right\\}$. Although the designed observer may not
be optimum in minimum $\mathcal{H}_{\infty}$ sense when using this alternative
method, we can guarantee the steady state convergence of the radii instead.

Algorithm 1 Fixed-Order Input & State Set-Valued Observer

1:Initialize: $M_{1}=\Sigma^{-1}$; $M_{2}=(C_{2}G_{2})^{\dagger}$;

2: $\Phi=I-G_{2}M_{2}C_{2}$;

3: Compute $\tilde{L}$ via Theorem 1;

4: $\Psi=I-\tilde{L}C_{2};$

5:
$\theta\triangleq\max_{i\in\\{1,2,\dots,N\\}}\|\Psi\Phi(A^{i}-G_{1}M_{1}C_{1})\|$;

6: $\hat{x}_{0|0}=\hat{x}_{0}=\text{centroid}(\hat{X}_{0})$;

7:
$\delta^{x}_{0}=\min\limits_{\delta}\\{\|x-\hat{x}_{0|0}\|\leq\delta,\forall
x\in\hat{X}_{0}\\}$;

8: $\hat{d}_{1,0}=M_{1}(z_{1,0}-C_{1}\hat{x}_{0|0}-D_{1}u_{0})$;

9:for $k=1$ to $K$ do

10:$\triangleright$ Estimation of $d_{2,k-1}$ and $d_{k-1}$

11:
$\hat{x}_{k|k-1}=\sum_{i=1}^{N}\lambda_{i,k}{A}^{i}\hat{x}_{k-1|k-1}+\sum_{i=1}^{N}\lambda_{i,k}{B}^{i}u_{k-1}$

12: $+G_{1}\hat{d}_{1,k-1};$

13:
$\hat{d}_{2,k-1}=M_{2}(z_{2,k}-{C}_{2}\hat{x}_{k|k-1}-\sum_{i=1}^{N}\lambda_{i,k}{D}^{i}_{2}u_{k});$

14: $\hat{d}_{k-1}=V_{1}\hat{d}_{1,k-1}+V_{2}\hat{d}_{2,k-1}$;

15:
$\delta^{d}_{k-1}=\delta_{k-1}^{x}\|V_{1}M_{1}C_{1}+V_{2}M_{2}C_{2}\hat{A}_{k}\|$

16:
$+\eta_{v}(\|(V_{2}M_{2}C_{2}G_{1}-V_{1})M_{1}T_{1}\|+\|V_{2}M_{2}T_{2}\|)$

17: $+\eta_{w}\|V_{2}M_{2}C_{2}\|;$

18:
$\hat{D}_{k-1}=\\{d\in\mathbb{R}^{l}:\|d-\hat{d}_{k-1}\|\leq\delta^{d}_{k-1}\\}$;

19:$\triangleright$ Time update

20: $\hat{x}^{\star}_{k|k}=\hat{x}_{k|k-1}+G_{2}\hat{d}_{2,k-1}$;

21:$\triangleright$ Measurement update

22:
$\hat{x}_{k|k}=\hat{x}^{\star}_{k|k}+\tilde{L}(z_{2,k}-C_{2}\hat{x}^{\star}_{k|k}-\sum_{i=1}^{N}\lambda_{i,k}D^{i}_{2}u_{k})$;

23:
$\delta^{x}_{k}=\delta^{x}_{0}\theta^{k}+\overline{\eta}\textstyle\sum_{i=1}^{k}\theta^{i-1}$;

24:
$\hat{X}_{k}=\\{x\in\mathbb{R}^{n}:\|x-\hat{x}_{k|k}\|\leq\delta^{x}_{k}\\}$;

25:$\triangleright$ Estimation of $d_{1,k}$

26:
$\hat{d}_{1,k}=M_{1}(z_{1,k}-C_{1}\hat{x}_{k|k}-\sum_{i=1}^{N}D^{i}_{1}u_{k})$;

27:end for

##  V Simulation Example

In this section, we consider a convex combination of two constituent linear
time-invariant strongly detectable subsystems that have been used in the
literature as a benchmark for some state and input filters (e.g., [12]):

$\displaystyle A^{1}=\begin{bmatrix}0.9&.5\\\
-0.3&1\end{bmatrix};A^{2}=\begin{bmatrix}0.85&.55\\\
-0.35&1\end{bmatrix};C=\begin{bmatrix}1&.2\\\ 1.1&1.9\end{bmatrix};$
$\displaystyle G=\begin{bmatrix}-0.02&0.04\\\
0.01&-0.05\end{bmatrix};H=\begin{bmatrix}1.1&2\\\
2.2&4\end{bmatrix};B^{1}=B^{2}=I_{2\times 2};D=0_{2\times 2}.$

The unknown inputs used in this example are as given in Figure 2, while the
initial state estimate and noise signals (drawn uniformly) have bounds
$\delta^{x}_{0}=0.5$, $\eta_{w}=0.02$ and $\eta_{v}=10^{-4}$. We also picked
uniformly random coefficients, $\lambda_{i,k}$, that satisfies
$0\leq\lambda_{i,k}\leq 1,\sum_{i=1}^{N}\lambda_{i,k}=1,\forall k$. Based on
the results of Theorem 1 and by solving the corresponding semi-definite
programming problem using YALMIP [15] and MOSEK [16] as the solver, we find
$S^{\star}=\begin{bmatrix}0.2745&0.1933\\\ 0.1933&0.4200\end{bmatrix}$,
$Y^{\star}=\begin{bmatrix}0.0010\\\ 0.1613\end{bmatrix}$ and the
$\mathcal{H}_{\infty}$-observer gain as
$\tilde{L}=S^{\star\,-1}Y^{\star}=\begin{bmatrix}-0.3946\\\
0.5656\end{bmatrix}$. Then, applying Algorithm 1, we summarized the set-valued
state and unknown input results in Figures 2 and 2. The radii are observed to
be convergent to steady state in Figure 2.

Figure 1: Actual states $x_{1}$, $x_{2}$ and their estimates, as well as
unknown inputs $d_{1}$ and $d_{2}$ and their estimates.

Figure 2: Actual estimation errors and radii of set-valued estimates of
states, $\|\tilde{x}_{k|k}\|$, $\delta^{x}_{k}$, and unknown inputs,
$\|\tilde{d}_{k}\|$, $\delta^{d}_{k}$.

##  VI Conclusion

We presented a fixed-order set-valued $\mathcal{H}_{\infty}$-observer for
linear parameter-varying bounded-error discrete-time dynamic systems, which
can be expressed as a convex combination of strongly detectable linear time-
invariant constituent systems. We provided sufficient conditions for the
optimality of the designed observer, which can be obtained from a semi-
definite programming problem with LMI constraints. We also showed that the
strong detectability of the constituent linear time-invariant systems is
necessary for the existence and stability of such an observer and for the
boundedness of the set-valued estimates. In addition, we came up with
sufficient structural conditions for the convergence of the radii of the set-
valued state and input estimates and derived the steady state radii. Finally,
we demonstrated the effectiveness of our proposed approach using an
illustrative example.

## References

  * [1] S. Gillijns and B. De Moor.  Unbiased minimum-variance input and state estimation for linear discrete-time systems.  Automatica, 43(1):111–116, January 2007.
  * [2] S. Gillijns and B. De Moor.  Unbiased minimum-variance input and state estimation for linear discrete-time systems with direct feedthrough.  Automatica, 43(5):934–937, 2007.
  * [3] S.Z. Yong, M. Zhu, and E. Frazzoli.  On strong detectability and simultaneous input and state estimation with a delay.  In IEEE Conference on Decision and Control (CDC), pages 468–475, 2015.
  * [4] S.Z. Yong, M. Zhu, and E. Frazzoli.  A unified filter for simultaneous input and state estimation of linear discrete-time stochastic systems.  Automatica, 63:321–329, 2016.
  * [5] P. Lu, E-J. van Kampen, C. C. de Visser, and Q. Chu.  Framework for state and unknown input estimation of linear time-varying systems.  Automatica, 73:145–154, 2016.
  * [6] S.Z Yong.  Simultaneous input and state set-valued observers with applications to attack-resilient estimation.  In 2018 Annual American Control Conference (ACC), pages 5167–5174. IEEE, 2018.
  * [7] M. Pajic, P. Tabuada, I. Lee, and G.J. Pappas.  Attack-resilient state estimation in the presence of noise.  In IEEE Conference on Decision and Control (CDC), pages 5827–5832, 2015.
  * [8] Y. Nakahira and Y. Mo.  Dynamic state estimation in the presence of compromised sensory data.  In IEEE Conference on Decision and Control (CDC), pages 5808–5813, 2015.
  * [9] S.Z. Yong, M.Q. Foo, and E. Frazzoli.  Robust and resilient estimation for cyber-physical systems under adversarial attacks.  In IEEE American Control Conference (ACC), pages 308–315, 2016\.
  * [10] F. Blanchini and M. Sznaier.  A convex optimization approach to synthesizing bounded complexity $\ell^{\infty}$ filters.  IEEE Transactions on Automatic Control, 57(1):216–221, 2012.
  * [11] M. Milanese and A. Vicino.  Optimal estimation theory for dynamic systems with set membership uncertainty: An overview.  Automatica, 27(6):997–1009, 1991.
  * [12] J. Chen and C.M. Lagoa.  Observer design for a class of switched systems.  In IEEE Conference on Decision and Control European Control Conference, pages 2945–2950, 2005.
  * [13] B.D.O. Anderson and J.B. Moore.  Detectability and stabilizability of time-varying discrete-time linear systems.  SIAM Journal on Control and Optimization, 19(1):20–32, 1981.
  * [14] M.A. Peters and P.A. Iglesias.  A spectral test for observability and reachability of time-varying systems.  SIAM Journal on Control Optimization, 37(5):1330–1345, August 1999\.
  * [15] J. Löfberg.  Yalmip : A toolbox for modeling and optimization in matlab.  In In Proceedings of the CACSD Conference, Taipei, Taiwan, 2004\.
  * [16] MOSEK ApS.  The MOSEK optimization toolbox for MATLAB manual. Version 8.1., 2017\.  URL: http://docs.mosek.com/8.1/toolbox/index.html.
  * [17] C.E de Souza, K.A Barbosa, and A.T. Neto.  Robust $\mathcal{H}_{\infty}$ filtering for discrete-time linear systems with uncertain time-varying parameters.  IEEE Transactions on Signal Processing, 54(6):2110–2118, 2006.

## Appendix: Proofs

###  VI-A Proof of Lemma 2

(11)-(15) and plugging $M_{1}=\Sigma^{-1}$ into (21) imply that

$\displaystyle\hat{d}_{1,k}$
$\displaystyle=M_{1}({C}_{1}\tilde{x}_{k|k}+\Sigma
d_{1,k}+\textstyle{\sum}_{i=1}^{N}\lambda_{i,k}{v}^{i}_{1,k}),$ (21)
$\displaystyle\hat{d}_{2,k-1}$
$\displaystyle=M_{2}(C_{2}(\textstyle{\sum}_{i=1}^{N}\lambda_{i,k-1}A^{i}\tilde{x}_{k-1|k-1}+G_{1}\tilde{d}_{1,k-1}$
(22) $\displaystyle\
+G_{2}d_{2,k-1}+\textstyle{\sum}_{i=1}^{N}\lambda_{i,k-1}w^{i}_{k-1})+\textstyle{\sum}_{i=1}^{N}\lambda_{i,k}v^{i}_{2,k}).$
$\displaystyle\tilde{d}_{1,k}$
$\displaystyle=d_{1,k}-\hat{d}_{1,k}=-M_{1}({C}_{1}\tilde{x}_{k|k}+\textstyle{\sum}_{i=1}^{N}\lambda_{i,k}{v}^{i}_{1,k}).$
(23)

(23) and setting $M_{2}=(C_{2}G_{2})^{\dagger}$ (Lemma 1) in (22), return

$\displaystyle\tilde{d}_{2,k-1}=-M_{2}(C_{2}\hat{A}_{k-1}\tilde{x}_{k-1|k-1}-C_{2}G_{1}M_{1}\textstyle{\sum}_{i=1}^{N}\lambda_{i,k-1}{v}^{i}_{1,k-1}$
$\displaystyle+\textstyle
C_{2}\textstyle{\sum}_{i=1}^{N}\lambda_{i,k-1}{w}^{i}_{k-1}+\textstyle{\sum}_{i=1}^{N}\lambda_{i,k}{v}^{i}_{2,k}).$
(24)

Defining $\tilde{x}^{\star}_{k|k}\triangleq x_{k}-\hat{x}^{\star}_{k|k}$, from
(1), (15) and (16) we obtain

$\displaystyle\tilde{x}^{\star}_{k|k}=\textstyle{\sum}_{i=1}^{N}\lambda_{i,k-1}({A}^{i}\tilde{x}_{k-1|k-1}+{w}^{i}_{k-1})+G_{1}\tilde{d}_{1,k-1}+\textstyle
G_{2}\tilde{d}_{2,k-1}$ (25)

In addition, from (11) and (19) and (23)-(25) we conclude:

$\displaystyle\tilde{x}_{k|k}$
$\displaystyle=(I-\tilde{L}C_{2})\tilde{x}^{\star}_{k|k}-\tilde{L}\textstyle{\sum}_{i=1}^{N}\lambda_{i,k}{v}^{i}_{2,k}.$
(26) $\displaystyle\tilde{x}^{\star}_{k|k}$
$\displaystyle=\overline{A}_{k-1}\tilde{x}_{k-1|k-1}-(I-G_{2}M_{2}C_{2})(G_{1}M_{1}$
(27)
$\displaystyle\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}(v^{i}_{1,k-1}-w^{i}_{k-1}))-G_{2}M_{2}\sum_{i=1}^{N}\lambda_{i,k}v^{i}_{2,k}.$

Now, we define  $\overline{w}_{k-1}\triangleq-
G_{2}M_{2}\textstyle\small{\sum_{i=1}^{N}}\lambda_{i,k}{v}^{i}_{2,k}-(I-G_{2}M_{2}C_{2})(G_{1}M_{1}\textstyle\small{\sum_{i=1}^{N}}\lambda_{i,k-1}({v}^{i}_{1,k-1}-{w}^{i}_{k-1})),\overline{v}_{k-1}\triangleq\textstyle\sum_{i=1}^{N}\lambda_{i,k}{v}^{i}_{2,k}$.
Then, (26)-(27) imply that

$\displaystyle\begin{array}[]{rl}\tilde{x}^{\star}_{k|k}&=\overline{A}_{k-1}\tilde{x}_{k-1|k-1}+\overline{w}_{k-1},\\\
\tilde{x}_{k|k}&=(I-\tilde{L}C_{2})\overline{A}_{k-1}\tilde{x}_{k-1|k-1}+(I-\tilde{L}C_{2})\overline{w}_{k-1}-\tilde{L}\overline{v}_{k-1}.\end{array}$
(30)

Now, consider the following linear time-varying system:

$\displaystyle
x_{k+1}=\overline{A}_{k}x_{k}+\overline{w}_{k},y_{k}=C_{2}x_{k}+\overline{v}_{k}.$
(31)

Systems (30) and (31) are equivalent from the viewpoint of estimation, since
the estimation error equations for both problems are the same, hence they both
have the same objective. Therefore, the pair $(\overline{A}_{k},C_{2})$ needs
to be uniformly detectable such that the observer is stable [13, Section 5]. ∎

###  VI-B Proof of Theorem 1

Starting from (30), we have

$\displaystyle\tilde{x}_{k|k}=(I-\tilde{L}C_{2})\overline{A}_{k-1}\tilde{x}_{k-1|k-1}+(I-\tilde{L}C_{2})\overline{w}_{k-1}-\tilde{L}\overline{v}_{k-1},$

from which we can define a system with $\tilde{x}_{k|k}$ as its state and
$\tilde{z}_{k|k}=\tilde{x}_{k|k}$ as the output:

$\displaystyle\begin{array}[]{rl}\tilde{x}_{k|k}&=(I-\tilde{L}C_{2})\overline{A}_{k-1}\tilde{x}_{k-1|k-1}+\begin{bmatrix}I-\tilde{L}C_{2}&-\tilde{L}\end{bmatrix}\begin{bmatrix}\overline{w}_{k-1}\\\
\overline{v}_{k-1}\end{bmatrix},\\\
\tilde{z}_{k|k}&=\tilde{x}_{k|k}.\end{array}$

By [17, Lemma 3], this system has an $\mathcal{H}_{\infty}$ performance
bounded by $\eta$, if there exists a symmetric positive definite matrix $P$
with rank $n$ such that:

$\displaystyle\begin{bmatrix}P&(I-\tilde{L}C_{2})\overline{{A}^{i}}P&\begin{bmatrix}I-\tilde{L}C_{2}&-\tilde{L}\end{bmatrix}&0\\\
*&P&0&P\\\ *&*&\eta I&0\\\ *&*&*&\eta I\end{bmatrix}\succ 0,\forall
i\in\left\\{1,2,\dots,N\right\\}.$ (32)

Notice that the referenced lemma requires the existence of a _bounded matrix
sequence_ , which in our case is a sequence of time-invariant matrices ($P$ is
the same for each $k$), that is obviously bounded. By plugging
$S={P}^{-1}\succ 0$ and applying some similarity transformations, we obtain

$\displaystyle\begin{bmatrix}0&S&0&0\\\ *&0&0&0\\\ *&*&I&0\\\
*&*&*&I\end{bmatrix}\begin{bmatrix}P&(I-\tilde{L}C_{2})\overline{{A}^{i}}P&\begin{bmatrix}I-\tilde{L}C_{2}&-\tilde{L}\end{bmatrix}&0\\\
*&P&0&P\\\ *&*&\eta I&0\\\ *&*&*&\eta I\end{bmatrix}\begin{bmatrix}0&S&0&0\\\
*&0&0&0\\\ *&*&I&0\\\ *&*&*&I\end{bmatrix}$
$\displaystyle=\begin{bmatrix}S&\overline{{A}^{i}}^{\top}(I-{C}_{2}^{\top}\tilde{L}^{\top})S&0&I\\\
*&S&\begin{bmatrix}I-\tilde{L}C_{2}&-\tilde{L}\end{bmatrix}&0\\\ *&*&I&0\\\
*&*&*&\eta I\end{bmatrix}\succ 0\ \forall i\in\left\\{1,2,\dots,N\right\\}.$

Setting $Y\triangleq S\tilde{L}$ completes the proof. ∎

###  VI-C Proof of Theorem 2

Suppose, for contradiction, that there exists an
$\mathcal{H}_{\infty}$-observer for system (1) with any sequence
$\\{\lambda_{i,k}\\}_{k=0}^{\infty}$ for all $i\in\\{1,2,\dots,N\\}$ that
satisfies $0\leq\lambda_{i,k}\leq 1,\sum_{i=1}^{N}\lambda_{i,k}=1,\forall k$,
but one of the constituent linear time-invariant systems (e.g.,
$({A}^{j},G,{C},H)$) is not strongly detectable. Since the
$\mathcal{H}_{\infty}$-observer exists for any sequence of $\lambda_{i,k}$,
particularly it exists when $\lambda_{j,k}=1$ and $\lambda_{ij,k}=0$, $\forall
i\neq j$ for all $k$. However, we know from [6] that strong detectability is
necessary for the stability of the linear time-invariant system
$({A}^{j},G,{C},H)$, which is a contradiction.∎

###  VI-D Proof of Theorem 3

To prove Theorem 3, we first find closed form expressions for the state and
input estimation errors in the following:

######  Lemma 3.

The state and input estimation errors are

$\displaystyle\tilde{x}_{k|k}$
$\displaystyle=(\textstyle\prod_{j=0}^{k-1}A_{e,k-j})\tilde{x}_{0|0}+\sum_{i=1}^{k}(\textstyle\prod_{j=0}^{i-2}A_{e,k-j})(\Psi\overline{w}_{k-i}-\tilde{L}\overline{v}_{k-i}),$
$\displaystyle\tilde{d}_{k-1}$
$\displaystyle=\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}(-V_{1}M_{1}C_{1}-V_{2}M_{2}C_{2}A_{e,i})\tilde{x}_{k-1|k-1}$
$\displaystyle+(-V_{1}M_{1}+V_{2}M_{2}C_{2}G_{1}M_{1})T_{1}\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}v^{i}_{k-1}$
$\displaystyle-
V_{2}M_{2}C_{2}\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}w^{i}_{k-1}-V_{2}M_{2}T_{2}\sum_{i=1}^{N}\lambda_{i,k}v^{i}_{k}.$

###### Proof.

From (30), we have

$\displaystyle\tilde{x}_{k|k}=\Psi\overline{A}_{k}\tilde{x}_{k-1|k-1}+\Psi\overline{w}_{k-1}-\tilde{L}\overline{v}_{k-1}.$
(33)

We use induction and (33) to obtain

$\displaystyle\begin{array}[]{rl}&\tilde{x}_{1|1}=\Psi\overline{A}_{1}\tilde{x}_{0|0}+\Psi\overline{w}_{0}-\tilde{L}\overline{v}_{0}=A_{e,1}\overline{x}_{0|0}+\Psi\overline{w}_{1-1}-\tilde{L}\overline{v}_{1-1}\\\
&=(\textstyle\prod_{j=0}^{1-1}A_{e,1-1})\tilde{x}_{0|0}+\textstyle\sum_{i=1}^{1}(\textstyle\prod_{j=0}^{i-2}A_{e,1-j})(\Psi\overline{w}_{k-1}-\tilde{L}\overline{v}_{1-i})\\\
&\tilde{x}_{k+1|k+1}=\Psi\overline{A}_{k+1}\tilde{x}_{k|k}+\Psi\overline{w}_{k}-\tilde{L}\overline{v}_{k}\\\
&=\Psi\overline{A}_{k+1}\big{[}(\textstyle\prod_{j=0}^{k-1}A_{e,k-j})\tilde{x}_{0|0}+\textstyle\sum_{i=1}^{k}(\textstyle\prod_{j=0}^{i-2}A_{e,k-j})(\Psi\overline{w}_{k-i}-\tilde{L}\overline{v}_{k-i})\big{]}\\\
&\quad+\Psi\overline{w}_{k}-\tilde{L}\overline{v}_{k}=(A_{e,k+1}A_{e,k}...A_{e,1})\tilde{x}_{0|0}+\Psi\overline{w}_{k}-\tilde{L}\overline{v}_{k}\\\
&\quad+\textstyle\sum_{i=1}^{k}(A_{e,k+1}A_{e,k}...A_{e,k-(i-2)})(\Psi\overline{w}_{k-i}-\tilde{L}\overline{v}_{k-i}))\\\
&=(\textstyle\prod_{j=0}^{k+1}A_{e,k+1-j})\tilde{x}_{0|0}+\textstyle\sum_{i=0}^{k}(\textstyle\prod_{j=0}^{i-2}A_{e,k-j})(\Psi\overline{w}_{k-i}-\tilde{L}\overline{v}_{k-i})\\\
&=(\textstyle\prod_{j=0}^{k+1}A_{e,k+1-j})\tilde{x}_{0|0}+\textstyle\sum_{i=1}^{k+1}(\textstyle\prod_{j=0}^{i-2}A_{e,k+1-j})(\Psi\overline{w}_{k+1-i}-\tilde{L}\overline{v}_{k+1-i}).\end{array}$

As for $\tilde{d}_{k-1}$, (23)-(VI-A) imply

$\displaystyle\begin{array}[]{ll}&\tilde{d}_{k-1}=V_{1}\tilde{d}_{1,k-1}+V_{2}\tilde{d}_{2,k-1}\\\
&=\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}(-V_{1}M_{1}C_{1}-V_{2}M_{2}C_{2}A_{e,i})\tilde{x}_{k-1|k-1}\\\
&+(-V_{1}M_{1}+V_{2}M_{2}C_{2}G_{1}M_{1})T_{1}\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}v^{i}_{k-1}\\\
&-V_{2}M_{2}C_{2}\textstyle\sum_{i=1}^{N}\lambda_{i,k-1}w^{i}_{k-1}-V_{2}M_{2}T_{2}\textstyle\sum_{i=1}^{N}\lambda_{i,k}v^{i}_{k}.\vspace{-0.4cm}\end{array}$
(38)

∎

Now, we are ready to prove Theorem 3. First, we define

$\displaystyle
B_{e,k}\triangleq\textstyle\prod_{j=0}^{k-1}A_{e,k-j},C^{i}_{e,k}\triangleq\textstyle\prod_{j=0}^{i-2}A_{e,k-j},\overline{t}_{k}\triangleq\Psi\overline{w}_{k}-\tilde{L}\overline{v}_{k}$
(39)

for $1\leq i\leq k$. Then, from Lemma 3, it follows that

$\displaystyle\begin{array}[]{rl}\|\tilde{x}_{k|k}\|&=\|B_{e,k}\tilde{x}_{0|0}+\textstyle\small{\sum_{i=1}^{k}}C^{i}_{e,k}\overline{t}_{k-i}\|\\\
&\leq\|B_{e,k}\|\|\tilde{x}_{0|0}\|+\|\textstyle\sum_{i=1}^{k}C^{i}_{e,k}\overline{t}_{k-i}\|.\end{array}$
(42)

Moreover, by similar reasoning,

$\displaystyle\|B_{e,k}\|=\|\textstyle\prod_{j=0}^{k-1}A_{e,k-j}\|\leq\textstyle\prod_{j=0}^{k-1}\|A_{e,k-j}\|$
$\displaystyle=\textstyle\prod_{j=0}^{k-1}\|\Psi\Phi\hat{A}_{k-j}\|=\textstyle\prod_{j=0}^{k-1}\|\Psi\Phi\textstyle\sum_{i=1}^{N}\lambda^{i}_{k-j}(A^{i}-G_{1}M_{1}C_{1})\|$
$\displaystyle=\textstyle\prod_{j=0}^{k-1}\|\textstyle\sum_{i=1}^{N}\lambda^{i}_{k-j}\Psi\Phi(A^{i}-G_{1}M_{1}C_{1})\|$
(43)
$\displaystyle\leq\textstyle\prod_{j=0}^{k-1}\textstyle\sum_{i=1}^{N}\lambda^{i}_{k-j}\|\Psi\Phi(A^{i}-G_{1}M_{1}C_{1})\|\leq\textstyle\prod_{j=0}^{k-1}\theta=\theta^{k},$
$\displaystyle\|$
$\displaystyle\textstyle\sum_{i=1}^{k}C^{i}_{e,k}\overline{t}_{k-i}\|\leq\textstyle\sum_{i=1}^{k}\|C^{i}_{e,k}\overline{t}_{k-i}\|\leq\textstyle\sum_{i=1}^{k}\|C^{i}_{e,k}\|\|\overline{t}_{k-i}\|,$
(44)
$\displaystyle\|C^{i}_{e,k}\|=\|\textstyle\prod_{j=0}^{i-2}A_{e,k-j}\|\leq\textstyle\prod_{j=0}^{i-2}\|A_{e,k-j}\|$
$\displaystyle=\textstyle\prod_{j=0}^{i-2}\|\textstyle\sum_{s=1}^{N}\lambda_{s,k-j}A_{e,s}\|\leq\textstyle\prod_{j=0}^{i-2}\theta\leq\theta^{i-1}.$
(45)

Furthermore, from the definition of $\overline{w}_{k}$ and (39) we have
$\overline{w}_{k-i}=-\Phi(G_{1}M_{1}\textstyle\sum_{s=1}^{N}\lambda_{s,k-i}v^{s}_{1,k-i}-\textstyle\sum_{s=1}^{N}\lambda_{s,k-i}w^{s}_{k-i})-G_{2}M_{2}\textstyle\sum_{s=1}^{N}\lambda_{s,k-i}v^{s}_{2,k-i},\|\overline{t}_{k-i}\|=\|\Psi\overline{w}_{k-i}-\tilde{L}\overline{v}_{k-i}\|=\|-\Psi\Phi
G_{1}M_{1}T_{1}\textstyle\sum_{s=1}^{N}\lambda_{s,k-i}v^{s}_{k-i}+\Psi\Phi\textstyle\sum_{s=1}^{N}\lambda_{s,k-i}w^{s}_{k-i}-\Psi
G_{2}M_{2}T_{2}\textstyle\sum_{s=1}^{N}\lambda_{s,k-i}v^{s}_{k-i}-\tilde{L}T_{2}\textstyle\sum_{s=1}^{N}v^{s}_{k-i}\|=\|\textstyle\sum_{s=1}^{N}\lambda_{s,k-i}(\Gamma
v^{s}_{k-i}+(\Psi\Phi)w^{s}_{k-i})\|\leq\overline{\eta}$,  from which, as well
as (42)-(45), we conclude that

$\displaystyle\|\tilde{x}_{k|k}\|\leq\|\tilde{x}_{0|0}\|\theta^{k}+\overline{\eta}\textstyle\sum_{i=1}^{k}\theta^{i-1}=\|\tilde{x}_{0|0}\|\theta^{k}+\overline{\eta}\frac{1-\theta^{k}}{1-\theta}\triangleq\delta^{x}_{k}.$
(46)

As for $\delta^{d}_{k-1}$, using Lemma 3 and (38), triangle inequality and the
facts that $0\leq\lambda_{i,k}\leq 1,\textstyle\sum_{i=1}^{N}\lambda_{i,k}=1$
and submultiplicativity of matrix norms, we obtain the result. ∎

###  VI-E Proof of Theorem 4

Notice that $0\leq\|A_{e,i}\|\leq\theta<1$ for all
$i\in\left\\{1,2,\dots,N\right\\}$ by assumption. So, $\theta^{k}$ in (46)
vanishes in steady state, which gives us the following steady state estimation
radius:
$\lim_{k\to\infty}\delta^{x}_{k}=\lim_{k\to\infty}\left(\|\tilde{x}_{0|0}\|\theta^{k}+\overline{\eta}\frac{1-\theta^{k}}{1-\theta}\right)=\frac{\overline{\eta}}{1-\theta}$.
Using this and starting from the expression for $\delta^{d}_{k-1}$ in Theorem
3, it converges to steady state, as follows:
$\lim_{k\to\infty}\delta^{d}_{k-1}=(\lim_{k\to\infty}\beta\delta^{x}_{k-1})+\|V_{2}M_{2}C_{2}\|\eta_{w}+(\|(V_{2}M_{2}C_{2}G_{1}-V_{1})M_{1}T_{1}\|+\|V_{2}M_{2}T_{2}\|)\eta_{v}=\frac{\overline{\eta}\beta}{1-\theta}+\eta_{w}\|V_{2}M_{2}C_{2}\|+\eta_{v}(\|V_{2}M_{2}T_{2}\|+\|R\|)$.∎


# Formation of Giant Clumps in High-$z$ Disc Galaxies by Compressive
Turbulence

Nir Mandelker,1 Omry Ginzburg,1 Avishai Dekel,1,2 Frederic Bournaud,3 Mark R.
Krumholz,4,5 Daniel Ceverino6,7 and Joel Primack8
1Centre for Astrophysics and Planetary Science, Racah Institute of Physics,
The Hebrew University, Jerusalem 91904, Israel
2SCIPP, University of California, Santa Cruz, CA 95064, USA
3AIM, CEA, CNRS, Universit´ e Paris-Saclay, Universit´ e Paris Diderot,
Sorbonne Paris Cit´ e, 91191 Gif-sur-Yvette, France
4Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
5Australian Research Council Centre of Excellence for All Sky Astrophysics in
3 Dimensions (ASTRO 3D), Australia
6Departamento de Fisica Teorica, Modulo 8, Facultad de Ciencias, Universidad
Autonoma de Madrid, 28049 Madrid, Spain
7CIAFF, Facultad de Ciencias, Universidad Autonoma de Madrid, 28049 Madrid,
Spain
8Department of Physics, University of California, Santa Cruz, CA, 95064, USA
E-mail<EMAIL_ADDRESS>
(Accepted –. Received –; in original form –)

###### Abstract

We address the formation of giant clumps in violently unstable gas-rich disc
galaxies at cosmic noon. While these are commonly thought to originate from
gravitational Toomre instability, cosmological simulations have indicated that
clumps form even in regions where the Toomre $Q$ parameter is well above
unity, which should be stable according to linear Toomre theory (Inoue et al.,
2016). Examining one of these cosmological simulations, we find that it
exhibits an excess in compressive modes of turbulence with converging motions.
The energy in converging motions within proto-clump regions is $\sim 70\%$ of
the total turbulent energy, compared to $\sim 17\%$ expected in equipartition.
When averaged over the whole disc, $\sim 32\%$ of the turbulent energy is in
converging motions, with a further $\sim 8\%$ in diverging motions. Thus, a
total of $\sim 40\%$ of the turbulent energy is in compressive modes, with the
rest in solenoidal modes, compared to the $(1/3)-(2/3)$ division expected in
equipartition. By contrast, we find that in an isolated-disc simulation with
similar properties, resembling high-$z$ star-forming galaxies, the energy in
the different turbulence modes are in equipartition, both in proto-clump
regions and over the whole disc. We conclude that the origin of the excessive
converging motions in proto-clump regions is external to the disc, and propose
several mechanisms that can induce them. This is an additional mechanism for
clump formation, complementary to and possibly preceding gravitational
instability.

###### keywords:

instabilities – methods: numerical – hydrodynamics – galaxies: evolution –
galaxies: formation

††pagerange: Formation of Giant Clumps in High-$z$ Disc Galaxies by
Compressive Turbulence–References††pubyear: 2018

##  1 Introduction

Roughly $60\%$ of $M_{\rm*}>10^{9}M_{\odot}$ star-forming galaxies (SFGs) at
$z\sim(1.5-3)$ are clumpy, with $\,{\rm kpc}$-scale clumps accounting for a
few percent of their mass and $\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10\%$ of their SFR (e.g. Guo et
al., 2015; Shibuya et al., 2016; Sattari et al., 2023). While a small fraction
of these clumps may be merging galaxies (Mandelker et al., 2014; Zanella et
al., 2019; Huško et al., 2023), the majority of them are understood to have
formed in situ through violent disc instabilities (VDI) (e.g. Noguchi, 1999;
Elmegreen et al., 2007; Bournaud et al., 2007; Genzel et al., 2008; Bournaud &
Elmegreen, 2009; Dekel et al., 2009; Ceverino et al., 2012; Mandelker et al.,
2014, 2017). VDI and clumpy star-formation are a robust feature of SFGs at
$z\sim 2$, in both observations and simulations. They are thought to play a
key role in mass flow to the galactic centre and the growth of bulges
(Elmegreen et al., 2008; Dekel et al., 2009; Ceverino et al., 2015; Zolotov et
al., 2015), the formation of dark matter cores (Ogiya & Nagai, 2022), the
chemical compositions of the thin and thick discs and the bulge (Clarke et
al., 2019; Debattista et al., 2023), globular cluster formation (Shapiro et
al., 2008; Mandelker et al., 2017), and turbulence driving in galactic discs
(Bournaud et al., 2009; Cacciato et al., 2012; Krumholz et al., 2018; Übler et
al., 2019; Ginzburg et al., 2022). The question of whether or not clumps
survive feedback and migrate towards the galactic center as bound structures
can constrain models of stellar feedback (Mandelker et al., 2017; Oklopčić et
al., 2017; Dekel et al., 2022, 2023a). Low-$z$ analogues of high-$z$ clumpy
discs have been observed in very gas rich systems (Fisher et al., 2017a, b;
Mehta et al., 2021), while analogues at cosmic dawn may be related to recent
observations of extremely bright galaxies with JWST (Dekel et al., 2023b).

The in situ formation of giant clumps is commonly attributed to Toomre
instability (Toomre, 1964; Dekel et al., 2009), where small fluctuations
collapse to clumps as their self-gravity ($\Sigma$) overcomes turbulent
pressure ($\sigma$) and rotational support ($\kappa$), namely once
$Q\propto\sigma\kappa/\Sigma<1$. Fragmentation into massive ($M_{\rm c}\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10^{8}{\rm M}_{\odot}$) long lived
clumps is commonly seen in simulations of isolated galaxy discs, especially
when the gas fraction is greater than $\sim 50\%$ (Fensch & Bournaud, 2021),
though Toomre-like instabilities manifest on multiple scales even when the gas
fraction is as low as $\sim 20\%$ (Renaud et al., 2021). Other studies of
idealized isolated galaxies suggest that while linear Toomre instability (or a
related spiral arm instability, Inoue & Yoshida, 2018, 2019) is responsible
for the initial fragmentation, this occurs on scales smaller than the observed
giant clumps which are themselves agglomerations of many small sub-clumps,
either physically bound or in projection (e.g Romeo & Agertz, 2014; Agertz et
al., 2015; Behrendt et al., 2015; Tamburello et al., 2015; Benincasa et al.,
2019). However, the Toomre $Q$ parameter is a linear concept applicable to
axisymmetric systems, while the violently unstable discs at high redshift are
far from being either linear or axisymmetric. Using the VELA suite of
cosmological simulations (Ceverino et al., 2014), which produce clumps with
properties consistent with observations (Mandelker et al., 2017; Guo et al.,
2018), Inoue et al. (2016) performed a detailed two-component (gas and stars)
analysis of the Toomre $Q$ parameter, following Romeo & Wiegert (2011). They
found, surprisingly, that massive clumps sometimes form in regions where
$Q>3$, which are expected to be stable to linear perturbations.

Several potential explanations for this have been proposed, such as the
importance of treating different gas phases separately in the $Q$ analysis
(Renaud et al., 2021), the multi-scale nature of $Q$ and of disc instabilities
more generally leading to small-scale fragmentation which then aggregate into
larger clumps (Romeo & Agertz, 2014; Agertz et al., 2015; Benincasa et al.,
2019), rapid decay of turbulence due to cooling and dissipation (Elmegreen,
2011), or various non-axisymmetric instabilities (Toomre, 1981; Inoue et al.,
2016; Inoue & Yoshida, 2018, 2019). However, as discussed in Inoue et al.
(2016), none of these seem sufficient to explain the simulation results.
Rather, they proposed that despite the results from isolated galaxy
simulations, VDI in realistic disc galaxies is not a simple linear
instability, but rather a highly perturbed configuration stimulated by a non-
linear process.

In this letter, we present a preliminary exploration of this possibility
focusing on the nature of turbulence in simulated discs. In particular, we
explore whether clump formation is associated with an enhancement in the ratio
of compressive to solenoidal (vortical) modes of turbulence. Previous work has
shown that such an enhancement can be associated with starbursts and the
formation of stellar clusters (Renaud et al., 2015a, b, 2022; Grisdale et al.,
2017). We explore this question using both cosmological and isolated galaxy
simulations, to examine whether the same mechanisms of clump formation are at
play in both cases. We introduce the simulations and our methods of clump
identification and turbulence decomposition in Section 2. In Section 3 we
present our results, and we discuss them and conclude in Section 4.

##  2 Method

###  2.1 Cosmological Simulation

We analyze one galaxy from the VELA suite of cosmological zoom-in simulations
(Ceverino et al., 2014), dubbed V07 (see table 1 of Mandelker et al., 2017).
The simulation utilised the ART code (Kravtsov et al., 1997; Kravtsov, 2003;
Ceverino & Klypin, 2009), which follows the evolution of a gravitating N-body
system and the Eulerian gas dynamics with an adaptive mesh refinement (AMR)
maximal resolution of $17.5-35$ physical $\,{\rm pc}$ at all times. The dark-
matter particle mass is $8.3\times 10^{4}{\rm M}_{\odot}$ and the minimal mass
of stellar particles is $10^{3}{\rm M}_{\odot}$. The code incorporates gas and
metal cooling, UV-background photoionization and self-shielding in dense gas,
stochastic star-formation, thermal feedback and metal enrichment from
supernovae and stellar winds (Ceverino et al., 2010, 2012) and feedback from
radiation pressure. Further details regarding the simulation method can be
found in Ceverino et al. (2014) and Mandelker et al. (2017) (models RadPRE and
RP, respectively).

The main VDI phase in V07 lasts from $z\sim 2.5-1$, following a dramatic phase
of “wet-compaction” to a star-forming “blue-nugget” (Zolotov et al., 2015).
During this time, the disc radius111See Mandelker et al. (2014, 2017) for
details on how $R_{\rm d}$ and other disc properties are defined in the
simulation. steadily increases from $R_{\rm d}\sim(9-20)\,{\rm kpc}$, its
thickness from $H_{\rm d}/R_{\rm d}\sim(0.1-0.2)$, and the disc crossing time
from $t_{\rm d}\sim R_{\rm d}/V_{\rm circ}\sim(30-60)\,{\rm Myr}$. The
baryonic disc mass increases from $M_{\rm d}\sim(2.7-6.5)\times 10^{10}{\rm
M}_{\odot}$, while the disc gas fraction decreases from $f_{\rm
g}\sim(0.28-0.13)$. The “cold” mass fraction in the disc, referring to cold
gas ($T<1.5\times 10^{4}\,{\rm K}$) and young stars ($<100\,{\rm Myr}$),
declines from $f_{\rm cold}\sim(0.4-0.14)$. V07 was analyzed in Inoue et al.
(2016) where it was shown that many giant clumps formed in regions with $Q>3$
and even as high as $Q\sim 10$. The properties of the clumps were analyzed in
detail in Mandelker et al. (2017) and their temporal evolution in Dekel et al.
(2022). Both of these suggest that the simulated clumps are consistent with
numerous observational constraints (see also Guo et al., 2018). As in Inoue et
al. (2016) and Dekel et al. (2022), we analyze the simulation with roughly 5
outputs per disc crossing time.

The left-hand panel of Fig. 1 shows a face-on projection of the gas surface
density at $z\sim 2.5$, near the onset of the VDI phase. The disc is fed by
two dense streams from the top and the bottom, which join the disc near its
radius at $R_{\rm d}\sim 12\,{\rm kpc}$, marked by a white circle. Clump
formation tends to occur in a thick ring near this radius (Dekel et al.,
2020).

###  2.2 Isolated Simulation

In addition to the cosmological simulation, we also analyse and idealised
isolated disc galaxy simulation that serves as a control case. This simulation
is similar to those described in Bournaud et al. (2014, e.g. Model G2), and in
Perret et al. (2014). The initial galaxy mass is $3.2\times 10^{10}{\rm
M}_{\odot}$ with an initial gas fraction of $65\%$, a disk half-mass radius of
$3.5\,{\rm kpc}$, and a disc-to-total ratio of $\sim 0.6$ inside the half-
stellar-mass radius. The simulation uses the RAMSES code (Teyssier, 2002), and
employs an adaptive-mesh grid with a maximal AMR resolution of $7.25\,{\rm
pc}$. The stellar feedback scheme includes supernovae feedback with thermal
energy injection and a sub-grid model for non-thermal processes (following
Teyssier et al., 2013), as well as radiation pressure and photo-ionization
from young massive stars (following Renaud et al., 2013).

Figure 1: Face on images of gas surface density in our cosmological simulation
at $z\sim 2.5$ (left), and our isolated galaxy simulation (right). Both images
are shortly after the onset of VDI. Dashed circles mark the disc radius,
$R_{\rm d}\sim 12\,{\rm kpc}$ and $3.5\,{\rm kpc}$ in the cosmological and
isolated simulations, respectively. The underlying white grid is meant to
guide the eye. The cosmological simulation image is integrated over $\pm
20\,{\rm kpc}$, and shows that the disc is fed by two streams from the top and
bottom, with most clump formation occurring in a thick ring near the radius
where the streams impact the disc. The isolated galaxy image is integrated
over $\pm 5\,{\rm kpc}$ and shows dense rings which are fragmenting into giant
clumps, as expected from Toomre instability.

The right-hand panel of Fig. 1 shows a face-on projection of the gas surface
density at $t\sim 140\,{\rm Myr}$ after the start of the simulation, shortly
after the disc begins to fragment into giant clumps. The gas mass fraction at
this time is $\sim 56\%$. The disc radius is $R_{\rm d}\sim 3.5\,{\rm kpc}$,
while clumps appear to form from the fragmentation of dense axisymmetric rings
at $R\sim(2-4)\,{\rm kpc}$.

###  2.3 Turbulence Decomposition

In general, any velocity field can be decomposed into a compressive part,
which has the same divergence as the initial field but no curl, and a
solenoidal part, which has the same curl as the initial field but no
divergence. For a turbulent (random) three-dimensional velocity field in
equipartition the solenoidal mode has twice as much power as the compressive
mode (Federrath et al., 2010), because the former consists of two dimensional
vortices while the latter consists of one dimensional shocks and expansion
waves. In general, the ratio of energies in each mode can deviate from this
value depending on the driver of the turbulence and the local dynamics. At a
given turbulent Mach number, the width of the density distribution increases
with increasing compressive to solenoidal ratio (Federrath et al., 2008),
which contributes to enhanced SFR in regions of excessive compressive modes
(e.g. Federrath & Klessen, 2012; Renaud et al., 2015a). Of particular
relevance to Toomre instability, if the velocity dispersion appearing in
$Q\propto\sigma\kappa/\Sigma$ consists of primarily compressive modes, in
particular converging flows, then a high value of $\sigma$ would not imply
gravitational stability but may actually trigger gravitational collapse. We
therefore wish to examine whether clump formation is associated with regions
of excessive compressive modes.

There are several methods to numerically evaluate the ratio of compressive to
solenoidal modes in a given region. Most studies of the turbulent ISM employ a
decomposition in Fourier space. However, such a decomposition can only be used
to estimate the total power in each mode in a given region and cannot easily
be used to study the spatial correlations of the mode ratios, because the
resulting compressive and solenoidal velocity fields are not locally
orthogonal in real-space (Brunt & Federrath, 2014). To get around this, one
can divide the full volume into sub-regions of a given size, perform a Fourier
analysis within each sub-region by erroneously assuming it to be periodic
rather than part of the full volume, and evaluate the local mode ratios (Vazza
et al., 2017). However, here we opt for a fully local approach in real-space.
A rigorous local approach is to decompose the shear tensor, the tensor of
first derivatives of the velocity field, around every point into three
components representing compression/expansion, solid body rotation, and shear
flow. We will discuss this approach and its relation to the Fourier
decomposition in an upcoming follow-up paper (Ginzburg et al., in
preparation). In this work, we instead take a simplified approach employed in
previous studies of the compressive to solenoidal ratio in galaxy-scale
simulations (Hopkins et al., 2013; Renaud et al., 2014, 2015a, 2015b; Grisdale
et al., 2017). We estimate the power in compressive modes at a given point by
$\sigma_{\rm c}^{2}\equiv|{\bf{\nabla}}\cdot{\bf{v}}|^{2}\Delta^{2}$, the
power in solenoidal modes by $\sigma_{\rm
s}^{2}\equiv|{\bf{\nabla}}\times{\bf{v}}|^{2}\Delta^{2}$, and the total
turbulent power as $\sigma_{\rm tot}^{2}\equiv\sigma_{\rm c}^{2}+\sigma_{\rm
s}^{2}$. Here, $\Delta$ is the length scale over which the divergence and curl
of the field are computed, though this cancels out when taking the ratio.
While $\sigma_{\rm c}$ and $\sigma_{\rm s}$ as defined above are not strictly
speaking ‘turbulence’, their ratio has been found to be a good proxy for the
ratio of turbulent modes. Following Renaud et al. (2014), we evaluate these on
a uniform grid with a grid spacing chosen to be at least 8 times larger than
the native simulation resolution, because the curl can artificially dissipate
on smaller scales (Federrath et al., 2011). For both cosmological and isolated
simulations, we deposit the gas onto a uniform grid with cell size
$\Delta=200\,{\rm pc}$ using a cloud-in-cell interpolation, and evaluate
derivatives using a three-point centered finite difference stencil. Using
$\Delta=400\,{\rm pc}$ or a two-point finite difference stencil does not
qualitatively change our results.

###  2.4 Clump Identification

We are primarily interested in correlating the ratio of compressive to
solenoidal power with the locations of protoclumps undergoing gravitational
collapse. We identify protoclumps following the same method employed in Inoue
et al. (2016) in their study of Toomre $Q$. We begin by identifying clumps in
the simulations using the method detailed in Mandelker et al. (2017). We
deposit the 3D density in a uniform grid with cell size $\Delta_{\rm cl}$, and
then smooth it with a spherical Gaussian with full-width-at-half-maximum $W$.
We then define $\delta\equiv(\rho-\rho_{{}_{\rm W}})/\rho$, with $\rho$ and
$\rho_{{}_{\rm W}}$ the unsmoothed and smoothed density fields, respectively.
We do this separately for the density of stars and of cold gas ($T<1.5\times
10^{4}\,{\rm K}$) plus young stars ($<100\,{\rm Myr}$), and take the larger of
the two $\delta$ values in each cell. We define clumps as connected regions
containing at least $8$ cells with $\delta>\delta_{\rm min}$. For the
cosmological simulation, we adopt $\Delta_{\rm cl}=70\,{\rm pc}$ ($2-4$ times
the maximal AMR resolution), $W=2.5\,{\rm kpc}$, and $\delta_{\rm min}=10$, as
in Mandelker et al. (2017) and Inoue et al. (2016). For the isolated
simulation, we use $\Delta_{\rm cl}=40\,{\rm pc}$ ($2$ times the maximal AMR
resolution), $W=0.5\,{\rm kpc}$, and $\delta_{\rm min}=15$. The different
values of $W$ and $\delta_{\rm min}$ used in the isolated simulation were
calibrated by eye based on face on images of the clumpy discs, and are due to
the more compact and gas rich nature of the isolated discs.

Clumps that contain at least 10 stellar particles are traced through time
based on their stellar particles. For each such clump at a given snapshot, we
search for all progenitor clumps in the preceding snapshot, defined as clumps
that contribute at least 25% of their stellar particles to the current clump.
If a given clump has more than one progenitor, we consider the most massive
one as the main progenitor. If a clump in snapshot $i$ has no progenitors in
snapshot $i-1$, we search the previous snapshots back to two disc crossing
times before snapshot $i$. If no progenitor is found in this period, snapshot
$i$ is declared the initial, formation time of the clump. We then trace all
clumps back in time one additional timestep to find the location of the
protoclump prior to collapse. We do so based on the center-of-mass velocity of
the clump at its initial time, in cylindrical coordinates in the disc frame.
Since we have $\sim 30$ snapshots per disc-orbital time in both simulations,
this extrapolation of the velocity field is reasonable. In what follows, we
limit our analysis to clumps with initial mass $M_{\rm c}>10^{7}{\rm
M}_{\odot}$, maximal mass over their lifetime $M_{\rm c,max}>10^{8}{\rm
M}_{\odot}$, formation site within $\pm 1.5\,{\rm kpc}$ of the disc midplane,
and survival time of at least one free-fall time, $t_{\rm
ff}\equiv\pi/2[R_{\rm c}^{3}/(GM_{\rm c})]^{1/2}$, with $M_{\rm c}$ the clump
baryonic mass and $R_{\rm c}$ the radius of a sphere with the same volume as
the clump.

##  3 Results

Figure 2: Fraction of total turbulent power in various modes in our
cosmological simulation. We show as a function of the cosmic expansion factor,
the mass-weighted average turbulent power within a cylinder of radius $R_{\rm
d}$ and half-thickness $1.5\,{\rm kpc}$. The blue (red) line is the fraction
of total power in solenoidal (compressive) modes. The cyan and green lines
further divide the compressive power into converging and expanding flows,
respectively. See text for details. On average, $\sim 60\%$ of the energy is
in solenoidal modes, close to the values of $\sim 2/3$ expected for
equipartition. Among the $\sim 40\%$ of the total energy in compressive modes,
most of this is in converging rather than expanding flows.  Figure 3: Clump
formation in the cosmological simulation. We show maps of gas surface density
(left) and the fraction of total turbulent energy in converging motions
(compressive modes with negative divergence, right). The maps have pixel sizes
of $200\,{\rm pc}$ and are integrated over $\pm 1.5\,{\rm kpc}$. We show the
evolution from $z\sim 2.49$ on the top to $z\sim 2.46$ on the bottom, with
roughly $\sim 9\,{\rm Myr}$ in between snapshots. We mark a forming clump from
its protoclump phase with a balck (white) circle in the left (right) columns.
The turbulence in the protoclump region is dominated by converging motions.

Figure 2 shows the fraction of total turbulent energy in various modes (see
Section 2.3) in our cosmological simulation, as a function of cosmic expansion
factor, $a$. For each mode, $\sigma^{2}$ represents the mass-weighted average
over a cylinder of radius $R_{\rm d}$ and half-thickness $1.5\,{\rm kpc}$,
where most clump formation occurs and where the gas densities are high. At
$z<3$ ($a>0.25$), $\sim 60\%$ of the turbulent energy is in solenoidal modes
(blue curve) while the remaining $\sim 40\%$ is in compressive modes (red
curve). This represents a slight ($\sim 10\%$) excess of compressive modes
compared to equipartition, where the expected fractions are $\sim 2/3$ and
$\sim 1/3$. At earlier times, there is a large excess of solenoidal modes,
presumably associated with the growth of angular momentum in the galaxy prior
to the formation of the disc. We further show the energy fraction in
converging (cyan) and expanding motions (green), defined as compressive modes
with negative and positive divergence, respectively. We find that $\sim 80\%$
of the energy in compressive modes is associated with converging motions,
compared to $\sim 20\%$ in expanding motions. This is consistent with the
compressive modes being dominated by regions of collapse rather than outflows.

In the isolated galaxy simulation (not shown), we find $\sim 70\%$ of the
energy in solenoidal modes and $\sim 30\%$ in compressive modes at all times,
very close to equipartition. Furthermore, most of the energy in compressive
modes is actually associated with expanding rather than converging motions.
The cosmological simulation thus exhibits an excess in compressive modes
associated with collapse compared to the isolated galaxies. This excess is
small when averaged over the entire disc, and the compressive modes never
dominate the global turbulence, as they may do following major-mergers
triggering galaxy-wide starbursts (Renaud et al., 2014). This is further
supported by the fact that the volume-weighted average values of $\sigma^{2}$
yield $\sim 68\%$ ($\sim 32\%$) in solenoidal (compressive) modes, consistent
with equipartition, with most of the compressive modes being associated with
expanding motions.

In Fig. 3, we examine the correlation between compressive modes, specifically
converging motions, and clump formation. We show a time sequence of our
cosmological simulation spanning $\sim 36\,{\rm Myr}$ from $z\sim(2.49-2.46)$.
We show the surface mass density of gas on the left, and the ratio of the
(mass-weighted) energy in converging motions to total turbulent energy along
the line-of-sight on the right. Each image has been integrated over $\pm
1.5\,{\rm kpc}$, as in Fig. 2, and has a pixel size of $200\,{\rm pc}$ (see
Section 2.3). In each panel, we highlight the position of a forming clump.
This clump is first identified at expansion factor $a=0.2872$ (second row),
and the position marked at $a=0.2866$ (top row) is the extrapolated protoclump
position (see Section 2.4). In the protoclump patch, converging motions
comprise $\sim 70\%$ of the total turbulent energy, compared to $\sim 17\%$
expected in equipartition where the energy in compressive modes is evenly
split between converging and expanding flows. This large excess of energy in
converging flows is maintained for the first $\sim 20\,{\rm Myr}$ of the clump
lifetime, after which it reduces to the equipartition value. It is worth
noting that after collapse, the clumps become dominated by rotation (Ceverino
et al., 2012) contributing to additional solenoidal modes inside the clump.
This same evolutionary sequence of clump formation was shown in Inoue et al.
(2016), Fig. 10, where it was shown that $Q\sim 5$ in the protoclump region,
and only reached $Q\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1$ after
$\sim 20\,{\rm Myr}$. We thus see that excessive compressive modes, and in
particular converging flows, can trigger clump formation in regions where
$Q\gg 1$.

Figure 4:  Distribution of the fraction of total turbulent energy in
converging motions in protoclumps (solid black and dotted blue) and their
antipodal regions (solid red and dotted magenta, see text) in the cosmological
(solid lines) and isolated galaxy simulations (dotted lines). The energies are
computed in a circular aperture of radius $500\,{\rm pc}$ centered on the
protoclumps or their antipodes, extending $\pm 1.5\,{\rm kpc}$ about the disc
midplane. In the isolated simulation, the protoclumps and the antipodes have
very similar distributions, sharply peaked at a ratio of $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}0.1$, slightly below the value
expected for equipartition. In the cosmological simulation, the protoclumps
have typical values of $\sim 0.6-0.7$ while the antipodes have typical values
of $\sim 0.3$. Excessive energy in converging flows is thus associated with
clump formation in the cosmological simulation, but not in the isolated galaxy
simulation.

One may worry that it is not the excessive converging motions that trigger
clump formation, but rather that the already collapsing protoclump induces
excessive converging flows. We address this, and systematically compare clump
formation in the cosmological and isolated galaxy simulations, in Fig. 4. We
show the distribution of the ratio of energy in converging motions to the
total turbulent energy in protoclump patches in the cosmological and isolated
galaxy simulation as solid black and dotted blue histograms, respectively. The
energies are mass-weighted averages within a circular aperture of radius
$500\,{\rm pc}$ around the center of the protoclump, and within $\pm 1.5\,{\rm
kpc}$ of the disc midplane as in Figs. 2 and 3. For each protoclump, we also
compute the energy ratio in the same aperture but centered on the antipodal
region, namely the region at the same galactocentric radius but on the
opposite side of the disc, $180^{\circ}$ from the protoclump. The distribution
of energy ratios in these antipodal regions are shown in solid red and dotted
magenta for the cosmological and isolated simulations, respectively.

In the isolated simulation, the ratio of energy in converging flows to total
turbulent energy is strongly peaked at $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}0.1$ in both the protoclumps and
their antipodes, slightly below the value expected for equipartition. This is
consistent with the global disc average where most of the compressive energy
is in expanding flows, presumably due to feedback. On the other hand, in the
cosmological simulations, the protoclump regions are strongly biased towards
large values of converging flows, typically $\sim 0.6-0.7$ of the total
turbulent energy. The converging-to-total energy ratios in the antipodes have
typical values of $\sim 0.3$, much lower than in the protoclumps though still
high compared to the isolated simulation and the equipartition value.

##  4 Discussion and conclusions

We have shown that violently unstable discs at $z\sim 2$ in cosmological
simulations exhibit an excess in compressive turbulence compared to the
equipartition value (Fig. 2). This excess is not seen in isolated galaxy
simulations. Furthermore, our results show that there is a qualitative
difference between the cosmological and isolated simulations in terms of the
correlation between clump formation and excessive compressive modes,
specifically converging motions with negative divergence. In cosmological
simulations, converging flows account for up to $\sim 70\%$ of the total
turbulent energy in protoclump regions, which can trigger clump formation even
in regions where $Q\gg 1$ (Fig. 3, compare to Fig. 10 in Inoue et al., 2016).
In isolated simulations, on the other hand, protoclump regions appear to be in
equipartition, with converging motions accounting for $\sim 1/6$ of the total
turbulent energy (Fig. 4). This suggests that the excessive energy in
converging flows seen within protoclumps in the cosmological simulation is of
an external origin; it is unlikely caused by the collapse of the clump itself,
or else we would have seen a similar signature in the isolated galaxy
simulation.

We conclude that the formation of giant clumps during VDI can be triggered by
converging flows associated with excessive compressive turbulence. This is a
non-linear process and is thus not directly related to the local value of
Toomre $Q$. The question is what causes such an excess. Since it is not seen
in isolated galaxy simulations, it must be external in origin. One possibility
is that it is induced by intense accretion from cold streams flowing along the
cosmic web. This is known to be an important aspect of VDI (Dekel et al.,
2009), and has been shown to be capable of inducing turbulence in discs at
levels comparable to both supernova feedback and gravitational instabilities
(Ginzburg et al., 2022), depending on the density distribution in the incoming
streams. As can be seen by comparing Figs. 1 and 3, the region with a large
excess of energy in converging motions where clumps tend to form correlates
with the region where inflowing streams join the disc. Another possibility is
that the excessive compressive modes are induced by an excess in fully
compressive tidal forces. This was shown to be the case in major-merger
simulations, where the tidal forces induced by the merger were fully
compressive across wide regions, leading to a large excess in compressive
turbulence which then triggered a starburst (Renaud et al., 2014, 2015a).
Similar tidally-induced clumps were found to form during major mergers in
cosmological simulations at $z\sim 8$ (Nakazato et al., 2024). The
cosmological simulation studied here does not undergo any major mergers, and
indeed the global ratio of compressive to solenoidal motions is only modestly
greater than equipartition. However, the chaotic nature of the galaxy
environment, including accretion of dense gas, minor mergers, pre-existing
clumps and other perturbations, may induce local tidal forces which are fully
compressive, leading to local excesses in converging flows and clump
formation. Finally, we cannot at this time rule out a third possibility, that
the difference is caused by the different feedback schemes implemented in the
cosmological and idealized simulations, rather than an external trigger.

Distinguishing these possibilities, and clarifying the origin of excessive
compressive modes in VDI galaxies more broadly, will be the focus of an
upcoming paper (Ginzburg et al., in preparation), where we examine a larger
sample of comsological and isolated galaxy simulations in more detail.

## Acknowledgements

We thank Frank van den Bosch, Christoph Federrath, Shigeki Inoue, and Florent
Renaud for very helpful and educational discussions. NM acknowledges support
from Israel Science Foundation (ISF) grant 3061/21. OG is supported by a
Milner fellowship and, along with AD, by ISF grant 861/20. DC is a Ramon-Cajal
Researcher and is supported by the Ministerio de Ciencia, Innovacion y
Universidades (MICIU/FEDER) under research grant PID2021-122603NB-C21. MRK
acknowledges support from the Australian Research Council through its Laureate
Fellowship program, award FL220100020, and from access to computational
resources supported by the Australian Government’s National Collaborative
Research Infrastructure Strategy and distributed through the National
Computational Merit Allocation Scheme, award jh2.

## References

  * Agertz et al. (2015) Agertz O., Romeo A. B., Grisdale K., 2015, MNRAS, 449, 2156
  * Behrendt et al. (2015) Behrendt M., Burkert A., Schartmann M., 2015, MNRAS, 448, 1007
  * Benincasa et al. (2019) Benincasa S. M., Wadsley J. W., Couchman H. M. P., Pettitt A. R., Tasker E. J., 2019, MNRAS, 486, 5022
  * Bournaud & Elmegreen (2009) Bournaud F., Elmegreen B. G., 2009, ApJ, 694, L158
  * Bournaud et al. (2007) Bournaud F., Elmegreen B. G., Elmegreen D. M., 2007, ApJ, 670, 237
  * Bournaud et al. (2009) Bournaud F., Elmegreen B. G., Martig M., 2009, ApJ, 707, L1
  * Bournaud et al. (2014) Bournaud F., et al., 2014, ApJ, 780, 57
  * Brunt & Federrath (2014) Brunt C. M., Federrath C., 2014, MNRAS, 442, 1451
  * Cacciato et al. (2012) Cacciato M., Dekel A., Genel S., 2012, MNRAS, 421, 818
  * Ceverino & Klypin (2009) Ceverino D., Klypin A., 2009, ApJ, 695, 292
  * Ceverino et al. (2010) Ceverino D., Dekel A., Bournaud F., 2010, MNRAS, 404, 2151
  * Ceverino et al. (2012) Ceverino D., Dekel A., Mandelker N., Bournaud F., Burkert A., Genzel R., Primack J., 2012, MNRAS, 420, 3490
  * Ceverino et al. (2014) Ceverino D., Klypin A., Klimek E. S., Trujillo-Gomez S., Churchill C. W., Primack J., Dekel A., 2014, MNRAS, 442, 1545
  * Ceverino et al. (2015) Ceverino D., Dekel A., Tweed D., Primack J., 2015, MNRAS, 447, 3291
  * Clarke et al. (2019) Clarke A. J., et al., 2019, MNRAS, 484, 3476
  * Debattista et al. (2023) Debattista V. P., et al., 2023, ApJ, 946, 118
  * Dekel et al. (2009) Dekel A., Sari R., Ceverino D., 2009, ApJ, 703, 785
  * Dekel et al. (2020) Dekel A., et al., 2020, MNRAS, 496, 5372
  * Dekel et al. (2022) Dekel A., Mandelker N., Bournaud F., Ceverino D., Guo Y., Primack J., 2022, MNRAS, 511, 316
  * Dekel et al. (2023a) Dekel A., Tziperman O., Sarkar K. C., Ginzburg O., Mandelker N., Ceverino D., Primack J., 2023a, MNRAS, 521, 4299
  * Dekel et al. (2023b) Dekel A., Sarkar K. C., Birnboim Y., Mandelker N., Li Z., 2023b, MNRAS, 523, 3201
  * Elmegreen (2011) Elmegreen B. G., 2011, ApJ, 737, 10
  * Elmegreen et al. (2007) Elmegreen D. M., Elmegreen B. G., Ravindranath S., Coe D. A., 2007, ApJ, 658, 763
  * Elmegreen et al. (2008) Elmegreen B. G., Bournaud F., Elmegreen D. M., 2008, ApJ, 688, 67
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Federrath et al. (2008) Federrath C., Klessen R. S., Schmidt W., 2008, ApJ, 688, L79
  * Federrath et al. (2010) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010, A&A, 512, A81
  * Federrath et al. (2011) Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, ApJ, 731, 62
  * Fensch & Bournaud (2021) Fensch J., Bournaud F., 2021, MNRAS, 505, 3579
  * Fisher et al. (2017a) Fisher D. B., et al., 2017a, MNRAS, 464, 491
  * Fisher et al. (2017b) Fisher D. B., et al., 2017b, ApJ, 839, L5
  * Genzel et al. (2008) Genzel R., et al., 2008, ApJ, 687, 59
  * Ginzburg et al. (2022) Ginzburg O., Dekel A., Mandelker N., Krumholz M. R., 2022, MNRAS, 513, 6177
  * Grisdale et al. (2017) Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093
  * Guo et al. (2015) Guo Y., et al., 2015, ApJ, 800, 39
  * Guo et al. (2018) Guo Y., et al., 2018, ApJ, 853, 108
  * Hopkins et al. (2013) Hopkins P. F., Narayanan D., Murray N., 2013, MNRAS, 432, 2647
  * Huško et al. (2023) Huško F., Lacey C. G., Baugh C. M., 2023, MNRAS, 518, 5323
  * Inoue & Yoshida (2018) Inoue S., Yoshida N., 2018, MNRAS, 474, 3466
  * Inoue & Yoshida (2019) Inoue S., Yoshida N., 2019, MNRAS, 488, 4400
  * Inoue et al. (2016) Inoue S., Dekel A., Mandelker, N. Ceverino D., Bournaud F., Primack J., 2016, MNRAS, 456, 2052
  * Kravtsov (2003) Kravtsov A. V., 2003, ApJ, 590, L1
  * Kravtsov et al. (1997) Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 73
  * Krumholz et al. (2018) Krumholz M. R., Burkhart B., Forbes J. C., Crocker R. M., 2018, MNRAS, 477, 2716
  * Mandelker et al. (2014) Mandelker N., Dekel A., Ceverino D., Tweed D., Moody C. E., Primack J., 2014, MNRAS, 443, 3675
  * Mandelker et al. (2017) Mandelker N., Dekel A., Ceverino D., DeGraf C., Guo Y., Primack J., 2017, MNRAS, 464, 635
  * Mehta et al. (2021) Mehta V., et al., 2021, ApJ, 912, 49
  * Nakazato et al. (2024) Nakazato Y., Ceverino D., Yoshida N., 2024, arXiv e-prints, p. arXiv:2402.08911
  * Noguchi (1999) Noguchi M., 1999, ApJ, 514, 77
  * Ogiya & Nagai (2022) Ogiya G., Nagai D., 2022, MNRAS, 514, 555
  * Oklopčić et al. (2017) Oklopčić A., Hopkins P. F., Feldmann R., Kereš D., Faucher-Giguère C.-A., Murray N., 2017, MNRAS, 465, 952
  * Perret et al. (2014) Perret V., Renaud F., Epinat B., Amram P., Bournaud F., Contini T., Teyssier R., Lambert J. C., 2014, A&A, 562, A1
  * Renaud et al. (2013) Renaud F., et al., 2013, MNRAS, 436, 1836
  * Renaud et al. (2014) Renaud F., Bournaud F., Kraljic K., Duc P. A., 2014, MNRAS, 442, L33
  * Renaud et al. (2015a) Renaud F., Bournaud F., Duc P.-A., 2015a, MNRAS, 446, 2038
  * Renaud et al. (2015b) Renaud F., et al., 2015b, MNRAS, 454, 3299
  * Renaud et al. (2021) Renaud F., Romeo A. B., Agertz O., 2021, MNRAS, 508, 352
  * Renaud et al. (2022) Renaud F., Segovia Otero Á., Agertz O., 2022, MNRAS, 516, 4922
  * Romeo & Agertz (2014) Romeo A. B., Agertz O., 2014, MNRAS, 442, 1230
  * Romeo & Wiegert (2011) Romeo A. B., Wiegert J., 2011, MNRAS, 416, 1191
  * Sattari et al. (2023) Sattari Z., et al., 2023, ApJ, 951, 147
  * Shapiro et al. (2008) Shapiro K. L., et al., 2008, ApJ, 682, 231
  * Shibuya et al. (2016) Shibuya T., Ouchi M., Kubo M., Harikane Y., 2016, ApJ, 821, 72
  * Tamburello et al. (2015) Tamburello V., Mayer L., Shen S., Wadsley J., 2015, MNRAS, 453, 2490
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Teyssier et al. (2013) Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068
  * Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
  * Toomre (1981) Toomre A., 1981, in Fall S. M., Lynden-Bell D., eds, Structure and Evolution of Normal Galaxies. pp 111–136
  * Übler et al. (2019) Übler H., et al., 2019, ApJ, 880, 48
  * Vazza et al. (2017) Vazza F., Jones T. W., Brüggen M., Brunetti G., Gheller C., Porter D., Ryu D., 2017, MNRAS, 464, 210
  * Zanella et al. (2019) Zanella A., et al., 2019, MNRAS, 489, 2792
  * Zolotov et al. (2015) Zolotov A., et al., 2015, MNRAS, 450, 2327