# Local and global canonical forms for differential-algebraic equations with
symmetries111Partially supported by the Research In Pairs program of
Mathematisches Forschungsinstitut Oberwolfach, whose hospitality is gratefully
acknowledged.

Peter Kunkel222 Mathematisches Institut, Universität Leipzig, Augustusplatz
10, D-04109 Leipzig, Fed. Rep. Germany<EMAIL_ADDRESS> Volker
Mehrmann333 Institut für Mathematik, TU Berlin, Str. des 17. Juni 136, D-10623
Berlin, Fed. Rep. Germany<EMAIL_ADDRESS>Partially supported by
the Deutsche Forschungsgemeinschaft through Project A2 of CRC 910 _Control of
self-organizing nonlinear systems: Theoretical methods and concepts of
application_

Abstract. Linear time-varying differential-algebraic equations with symmetries
are studied. The structures that we address are self-adjoint and skew-adjoint
systems. Local and global canonical forms under congruence are presented and
used to classify the geometric properties of the flow associated with the
differential equation as symplectic or generalized orthogonal flow. As
applications, the results are applied to the analysis of dissipative
Hamiltonian systems arising from circuit simulation and incompressible flow.

Keywords. Differential-algebraic equation, self-adjoint system, skew-adjoint
system, dissipative Hamiltonian system, canonical form under congruence.

AMS(MOS) subject classification. 37J05, 65L80, 65L05, 65P10, 49K15

Dedicated to Alfi Quarteroni on the occasion of his 70th birthday.

##  1 Introduction

We study regular linear variable coefficient differential-algebraic equations
(DAEs)

$E(t)\dot{x}=A(t)x+f(t),\quad\begin{array}[]{l}E,A\in
C({\mathbb{I}},{\mathbb{R}}^{n,n}),\ f\in
C({\mathbb{I}},{\mathbb{R}}^{n})\mbox{ sufficiently smooth},\end{array}$ (1)

which additionally possess certain symmetries, in particular self-adjoint and
skew-adjoint structures. Here ${\mathbb{I}}\subseteq{\mathbb{R}}$ is a compact
non-trivial time interval and $C^{k}({\mathbb{I}},{\mathbb{R}}^{n,n})$ with
$k\in{\mathbb{N}}_{0}\cup\\{\infty\\}$ denotes the set of $k$ times
continuously differentiable functions from $\mathbb{I}$ into the set of real
$n\times n$ matrices ${\mathbb{R}}^{n,n}$. In the case $k=0$ we drop the
superscript. A function is said to be sufficiently smooth if $k$ is
sufficiently large such that all needed derivatives exist.

The discussed classes of DAEs with symmetries are defined as follows.

###### Definition 1.1

The DAE (1) and its associated pair $(E,A)$ of matrix functions are called
self-adjoint if

$E^{T}=-E,\quad A^{T}=A+\dot{E}$ (2)

as equality of functions.

###### Definition 1.2

The DAE (1) and its associated pair $(E,A)$ of matrix functions are called
skew-adjoint if

$E^{T}=E,\quad A^{T}=-A-\dot{E}$ (3)

as equality of functions.

Systems with these type of symmetries arise in the modeling of physical
systems, e.g., by a (generalized) Hamiltonian formalism or in optimal control
problems leading to a self-adjoint structure or by so-called dissipative
Hamiltonian systems leading to a skew-adjoint structure.

Our main motivation to study differential-algebraic equations with self-
adjoint and skew-adjoint pairs arises from multi-physics, multi-scale models
that are coupling different physical domains that may include mechanical,
mechatronic, fluidic, thermic, hydraulic, pneumatic, elastic, plastic, or
electric components, see e.g. [1, 2, 9, 17, 35, 36].

An important class of problems where multi-physics and multi-scale modeling
arises is the human cariovascular system, where model hierarchies of detailed
models for the blood flow in large vessels, modeled via the Navier-Stokes
equations, and reduced or surrogate models for the capillary vessels, modeled
via electrical network equations, are coupled together to improve
computational efficiency, while at the same time achieving a desired
simulation accuracy, see [11, 12, 27, 28, 30]. Due to the physical background,
after space discretization and linearization along a solution, as well as
ignoring dissipation terms, all the components arising in this application can
be expressed as DAE systems with symmetries.

To deal with general multi-physics and multi-scale coupling in recent years
the framework of port-Hamiltonian systems has become an important modeling
paradigm [3, 13, 18, 25, 26, 33, 34] that encodes underlying physical
principles directly into the algebraic structure of the coefficient matrices
and in the geometric structure associated with the flow of the dynamical
system. This leads to a remarkably flexible modeling approach, which has also
been extended to include algebraic constraints, so that the resulting model is
a port-Hamiltonian differential-algebraic equation (pHDAE), [3, 25, 39]. Such
systems allow for automated modeling in a modularized network based fashion,
and they are ideal for building model hierarchies of very fine models for
numerical simulation and reduced or surrogate models for control and
optimization. This makes them particularly suited also for large networks,
such as power, gas, or district heating networks where such model hierarchies
are used to adapt the simulation and optimization techniques to user needs,
[7, 16, 25].

Since in this paper we will mainly discuss linear time-varying DAE systems
with symmetries, we introduce the structure of pHDAEs for this case as in [3],
see [25] for the more general nonlinear framework.

Linear time-varying _pHDAE systems with quadratic Hamiltonian_ have the form

$\displaystyle E(t)\dot{x}+E(t)K(t)x$
$\displaystyle=(J(t)-R(t))x+(G(t)-P(t))u,$ (4a) $\displaystyle y$
$\displaystyle=(G(t)+P(t))^{T}x+(S(t)-N(t))u,$ (4b)

with state $x$, input $u$, output $y$ and coefficients $E\in
C^{1}(\mathbb{I},\mathbb{R}^{n,n})$, $J,R,K\in
C(\mathbb{I},\mathbb{R}^{n,n})$, $G,P\in C(\mathbb{I},\mathbb{R}^{n,m})$,
$S,N\in C(\mathbb{I},\mathbb{R}^{m,m})$, $S=S^{T},N=-N^{T}$. As energy
function one has the quadratic _Hamiltonian_

$\mathcal{H}\colon\mathbb{I}\times\mathbb{R}^{n}\to\mathbb{R},\qquad(t,x)\mapsto\tfrac{1}{2}x^{H}E(t)x,$
(5)

and the pair of coefficients $(E,J-EK)$ is skew-adjoint, while the matrix
function associated with dissipation of energy

$\mathcal{W}=\begin{bmatrix}R&P\\\ P^{T}&S\end{bmatrix}$

is pointwise positive semidefinite. Furthermore, typically one also has that
$E$ is pointwise positive semidefinite as well. If the system does not have an
output equation and the input is considered as an inhomogeneity then this is
called a dissipative Hamiltonian DAE (dHDAE).

The underlying skew-adjoint structure arises if the dissipation term $R$ is
neglected, i.e., if $\mathcal{W}=0$ in (4). Typically the problems with
dissipation can be considered as a perturbed symmetry structure and a
dissipative term can be treated separately as a by-product in simulation
methods.

To illustrate applications for dHDAEs and DAEs with symmetries consider the
following simple examples, for further applications see [3, 25, 39, 40].

###### Example 1.3

Consider the pHDAE formulation of an electrical circuit from [25].

$+$$-$$E_{G}$$R_{G}$$I_{G}$$L$$I$$R_{L}$$R_{R}$$I_{R}$$C_{1}$$V_{1}$$I_{1}$$C_{2}$$V_{2}$$I_{2}$
Figure 1: circuit example

Denoting by $V_{i}$ the voltages and by $I_{i}$ the currents, where $L>0$
models an inductor, $C_{1},C_{2}>0$ capacitors, $R_{G},R_{L},R_{R}>0$
resistances, and $E_{G}$ a controlled voltage source, one obtains a pHDAE

$\displaystyle E\dot{x}$ $\displaystyle=(J-R)x+Gu,$ (6a) $\displaystyle y$
$\displaystyle=G^{T}x,$ (6b)

with $x=\begin{bmatrix}I&V_{1}&V_{2}&I_{G}&I_{R}\end{bmatrix}^{T}$, $u=E_{G}$,
$y=I_{G}$, $E=\mathop{\rm diag}\nolimits(L,C_{1},C_{2},0,0)$,

$G=\begin{bmatrix}0\\\ 0\\\ 0\\\ 1\\\ 0\end{bmatrix},\quad
J=\begin{bmatrix}0&-1&1&0&0\\\ 1&0&0&-1&0\\\ -1&0&0&0&-1\\\ 0&1&0&0&0\\\
0&0&1&0&0\end{bmatrix},\quad R=\begin{bmatrix}R_{L}&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&0\\\ 0&0&0&R_{G}&0\\\ 0&0&0&0&R_{R}\end{bmatrix}.$

The quadratic Hamiltonian, describing the energy stored in the inductor and
the two capacitors, is given by

$\mathcal{H}(I,V_{1},V_{2})=\frac{1}{2}LI^{2}+\frac{1}{2}C_{1}V_{1}^{2}+\frac{1}{2}C_{2}V_{2}^{2}.$

If the generator is shut down (i.e. $E_{G}=0$), then the system approaches an
equilibrium solution for which $\vphantom{\tilde{E}}{\textstyle{d\over
dt}}\mathcal{H}(x)=0$, so that $I=I_{G}=I_{R}=0$ and thus $x=0$ in the
equilibrium. Without the damping by the resistances (i.e. setting
$R_{L}=R_{G}=R_{R}=0$) this is a skew-adjoint DAE.

###### Example 1.4

A classical example of a partial differential equation which, after proper
space discretization leads to a pHDAE, see e.g. [10, 29], are the
incompressible or nearly incompressible Navier-Stokes equations describing the
flow of a Newtonian fluid in a domain $\Omega$,

$\displaystyle\partial_{t}{v}$ $\displaystyle=\nu\Delta
v-(v\cdot\nabla)v-\nabla p+f\qquad\text{in }\Omega\times\mathbb{T},$
$\displaystyle 0$ $\displaystyle=\nabla^{T}v\qquad\text{in
}\Omega\times\mathbb{T},$

together with suitable initial and boundary conditions, see e.g. [38]. When
one linearizes around a prescribed vector field $v_{\infty}$, then one obtains
the _linearized Navier-Stokes equations_ , and if $v_{\infty}$ is constant in
space and the term $(v_{\infty}\cdot\nabla)v$ is neglected then one obtains
the _Stokes equation_.

Performing a finite element discretization in space, see for instance [23], a
Galerkin projection leads to a dHDAE of the form

$\begin{bmatrix}M&0\\\ 0&0\end{bmatrix}\begin{bmatrix}\dot{v}\\\
\dot{p}\end{bmatrix}=\begin{bmatrix}A_{S}(t)-A_{H}(t)&-B\\\
B^{T}&-C\end{bmatrix}\begin{bmatrix}v\\\ p\end{bmatrix}+\begin{bmatrix}f(t)\\\
0\end{bmatrix},$ (7)

where $M=M^{T}>0$ is the mass matrix, $A_{S}=-A_{S}^{T}$, $A_{H}=A_{H}^{T}\geq
0$ are skew-symmetric and symmetric part of the discretized and linearized
convection-diffusion operator, $B$ is the discretized gradient operator,
$B^{T}$ is the discretized divergence operator, which we assume to be
normalized so that it is of full row rank, and $C=C^{T}>0$ is a stabilization
term of small norm that is needed for some finite element spaces, see e.g. [4,
29, 31, 32]. Here $v$ and $p$ denote the discretized velocity and pressure,
respectively, and $f$ is a forcing term. Without the damping (i.e. for
$A_{H}=0$ and $C=0$) we again have a skew-adjoint DAE.

The class of problems with self-adjoint structure arises most prominently in
the context of constrained generalized Hamiltonian systems and in optimal
control problems, where the operators associated with the optimality
conditions have this structure.

###### Example 1.5

In [21] the optimality conditions were derived for the linear-quadratic
optimal control problem of minimizing a cost functional

${\mathcal{J}}(x,u)=\frac{1}{2}x(t_{f})^{T}M_{f}x(t_{f})+\frac{1}{2}\int_{t_{0}}^{t_{f}}\left(x^{T}Wx+x^{T}Su+u^{T}S^{T}x+u^{T}Ru\right)\,dt,$
(8)

subject to the DAE constraint

$E\dot{x}=Ax+Bu+f,\quad x(t_{0})={x_{0}}\in{\mathbb{R}}^{n},$ (9)

with $E,A,W\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$, $B\in
C({\mathbb{I}},{\mathbb{R}}^{n,m})$, $S\in
C({\mathbb{I}},{\mathbb{R}}^{n,m})$, $R\in
C({\mathbb{I}},{\mathbb{R}}^{m,m})$, $f\in C({\mathbb{I}},{\mathbb{R}}^{n})$
and $M_{e}\in{\mathbb{R}}^{n,n}$, where $R=R^{T}$, $W=W^{T}$ and
$M_{f}=M_{f}^{T}$.

After some appropriate reformulation (via some index reduction process) and
under some smoothness conditions, the optimality condition is given by a
boundary value problem

$\left[\begin{array}[]{ccc}0&E&0\\\ -E^{T}&0&0\\\
0&0&0\end{array}\right]\frac{d}{dt}\left[\begin{array}[]{c}\lambda\\\ x\\\
u\end{array}\right]=\left[\begin{array}[]{ccc}0&A&B\\\
A^{T}+\frac{d}{dt}E^{T}&W&S\\\
B^{T}&S^{T}&R\end{array}\right]\left[\begin{array}[]{c}\lambda\\\ x\\\
u\end{array}\right]+\left[\begin{array}[]{c}f\\\ 0\\\ 0\end{array}\right],$

with boundary conditions $x(t_{0})={x_{0}}$,
$E(t_{f})^{T}\lambda(t_{f})-M_{f}x(t_{f})=0$. The associated pair of
coefficient functions obviously is a a self-adjoint pair, see [22].

###### Example 1.6

Linear multibody systems with linear holonomic constraints, see [14, 24], take
the form

$\displaystyle M\dot{p}$ $\displaystyle=$ $\displaystyle-Wq-G^{T}\lambda,$
$\displaystyle\dot{q}$ $\displaystyle=$ $\displaystyle p,$ $\displaystyle 0$
$\displaystyle=$ $\displaystyle Gq,$

where $p,q\in C^{1}(\mathbb{I},\mathbb{R}^{n})$, $W,M\in\mathbb{R}^{n,n}$ with
$W=W^{T}$, $M=M^{T}$, and $G\in\mathbb{R}^{m,n}$. If the mass matrix is
positive definite, i.e. $M>0$, then we can multiply the second equation by
$-M$ and the constraint by $-1$ to obtain, after switching the first and
second equation,

$\left[\begin{array}[]{ccc}0&M&0\\\ -M&0&0\\\
0&0&0\end{array}\right]\left[\begin{array}[]{c}\dot{q}\\\ \dot{p}\\\
\dot{\lambda}\end{array}\right]=\left[\begin{array}[]{ccc}-W&0&-G^{T}\\\
0&-M&0\\\ -G&0&0\end{array}\right]\left[\begin{array}[]{c}q\\\ p\\\
\lambda\end{array}\right],$

which has the structure of a self-adjoint DAE.

If $W>0$ then we can also multiply the second equation of the constrained
Hamiltonian system by $W$ to obtain

$\left[\begin{array}[]{ccc}W&0&0\\\ 0&M&0\\\
0&0&0\end{array}\right]\left[\begin{array}[]{c}\dot{q}\\\ \dot{p}\\\
\dot{\lambda}\end{array}\right]=\left[\begin{array}[]{ccc}0&W&0\\\
-W&0&-G^{T}\\\ 0&G&0\end{array}\right]\left[\begin{array}[]{c}q\\\ p\\\
\lambda\end{array}\right]$

which now has the structure of a skew-adjoint DAE.

###### Remark 1.7

It should be noted that if in a self-adjoint system $E\dot{x}=Ax+f$, both $E$
and $A$ are constant in time and invertible, then by multiplication with
$E^{-1}$ and a change of variables $z=Ax$ we get a system
$A^{-1}\dot{z}=E^{-1}z+E^{-1}f$ that is skew-adjoint. A similar construction
has also been proposed for dissipative Hamiltonian systems in [8].

Having illustrated the importance of DAEs with symmetries, in this paper we
present canonical forms for DAEs (1) with self-adjoint and skew-adjoint
structure and show the consequences for the resulting flows.

The paper is organized as follows. In Section 2 we recall some results for
general DAEs. In Section 3 we discuss canonical forms for self-adjoint and in
Section 4 for skew-adjoint DAEs.

##  2 Preliminaries

Linear time-varying DAE systems have been extensively discussed, see [20] for
a detailed analytical and numerical treatment. In this section we recall some
basic concepts from the general theory of DAEs. Our first concept is that of
_regularity_ for DAEs that is concerned with the existence of solutions at
least for sufficiently smooth inhomogeneity (surjectivity) and uniqueness of
the solution in the cases where the initial condition $x(t_{0})=x_{0}$ allows
for a solution (injectivity).

###### Definition 2.1

The pair $(E,A)$ and the corresponding DAE (1) are called _regular_ if

  1. 1.

the DAE (1) is solvable for every sufficiently smooth $f$,

  2. 2.

the solution is unique for every $t_{0}\in{\mathbb{I}}$ and every
$x_{0}\in{\mathbb{R}}^{n}$ allowing for a solution of the DAE with
$x(t_{0})=x_{0}$,

  3. 3.

the solution depends smoothly on $f$, $t_{0}$, and $x_{0}$.

The most important technique for the analysis of general linear DAEs is the
construction of suitable local and global canonical forms under global
equivalence. Since we will refer to these results and some techniques from
their derivation, we include here the necessary material from the general
(square) case.

We start with (global) equivalence which refers to time-dependent scaling of
the DAE and changes of basis.

###### Definition 2.2

Two pairs $(E_{i},A_{i})$, $E_{i},A_{i}\in
C({\mathbb{I}},{\mathbb{R}}^{n,n})$, $i=1,2$, of matrix functions are called
$($globally$)$ equivalent if there exist pointwise nonsingular matrix
functions $P\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$ and $Q\in
C^{1}({\mathbb{I}},{\mathbb{R}}^{n,n})$ such that

$E_{2}=PE_{1}Q,\quad A_{2}=PA_{1}Q-PE_{1}\dot{Q}$ (10)

as equality of functions. We write $(E_{1},A_{1})\sim(E_{2},A_{2})$.

It is easy to see that the relation defined in Definition 2.2 indeed is an
equivalence relation, [20].

The derivation of canonical forms then relies on the following property of
matrix functions on intervals, see [6, 20].

###### Theorem 2.3

Let $E\in C^{k}({\mathbb{I}},{\mathbb{R}}^{m,n})$,
$k\in{\mathbb{N}}_{0}\cup\\{\infty\\}$, with $\mathop{\rm rank}\nolimits
E(t)=r$ for all $t\in{\mathbb{I}}$. Then there exist pointwise orthogonal
functions $U\in C^{k}({\mathbb{I}},{\mathbb{R}}^{m,m})$ and $V\in
C^{k}({\mathbb{I}},{\mathbb{R}}^{n,n})$ such that

$U^{T}EV=\begin{bmatrix}\Sigma&0\\\ 0&0\end{bmatrix}$ (11)

with pointwise nonsingular $\Sigma\in C^{k}({\mathbb{I}},{\mathbb{R}}^{r,r})$.

We then have the following result on a local canonical form, i. e., a
canonical form that requires the restriction to certain subintervals, see [20,
Section 3.1,Section 3.3]. In the following, non-specified blocks of matrices
or matrix functions are denoted by $*$.

###### Theorem 2.4

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth. Then there exist pairwise disjoint open intervals
${\mathbb{I}}_{j}$, $j\in{\mathbb{N}}$, with

$\overline{\bigcup_{j\in{\mathbb{N}}}{\mathbb{I}}_{j}}={\mathbb{I}}$ (12)

such that on every ${\mathbb{I}}_{j}$ one has

$(E,A)\sim\left(\left[\begin{array}[]{cc}I_{d}&W\\\
0&G\end{array}\right],\left[\begin{array}[]{cc}0&0\\\
0&I_{a}\end{array}\right]\right),$ (13)

where $G$ is structurally nilpotent according to

$G=\left[\begin{array}[]{cccc}0&*&\cdots&*\\\ &\ddots&\ddots&\vdots\\\
&&\ddots&*\\\ &&&0\end{array}\right].$ (14)

Furthermore, the size $d$ of the differential part and the size $a$ of the
algebraic part are the same for every interval.

A global canonical form, i.e., a canonical form that does not require the
restriction to certain subintervals, was given in [5]. We state this result
here in a version for real-valued problems omitting the last step of the proof
which would require complex-valued transformations. We include the proof for
later reference.

###### Theorem 2.5

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth. Then we have

$(E,A)\sim\left(\left[\begin{array}[]{cc}I_{d}&E_{12}\\\
0&E_{22}\end{array}\right],\left[\begin{array}[]{cc}0&A_{12}\\\
0&A_{22}\end{array}\right]\right),$ (15)

where

$E_{22}(t)\dot{x}_{2}=A_{22}(t)x_{2}+f_{2}(t)$ (16)

is uniquely solvable for every sufficiently smooth $f_{2}$ without specifying
initial conditions.

Proof. If the homogeneous equation

$E(t)\dot{x}=A(t)x$

has only the trivial solution, then the first block in (15) is missing (i. e.,
$d=0$) and the claim holds trivially by assumption. In any case, the solution
space is finite dimensional, since otherwise we could not select a unique
solution by prescribing initial conditions.

Let therefore $d\neq 0$ and let $\\{\phi_{1},\ldots,\phi_{d}\\}$ be a basis of
the solution space. Setting $\Phi=[\>\phi_{1}\>\>\cdots\>\>\phi_{d}\>]$, we
have

$\mathop{\rm rank}\nolimits\Phi(t)=d$ for all $t\in{\mathbb{I}}$,

since, if we had $\mathop{\rm rank}\nolimits\Phi(t)<d$ for some
$t_{0}\in{\mathbb{I}}$, then there would exist coefficients
$\alpha_{1},\ldots,\alpha_{d}\in{\mathbb{R}}$, not all being zero, with

$\alpha_{1}\phi_{1}(t_{0})+\cdots+\alpha_{d}\phi_{d}(t_{0})=0$

and $\alpha_{1}\phi_{1}+\cdots+\alpha_{d}\phi_{d}$ would be a nontrivial
solution of the homogeneous initial value problem.

Hence, by Theorem 2.3 there exists a smooth, pointwise nonsingular matrix
function $U$ with

$U^{H}\Phi=\left[\begin{array}[]{c}I_{d}\\\ 0\end{array}\right].$

Defining

$\Phi^{\prime}=U\left[\begin{array}[]{c}0\\\ I_{a}\end{array}\right]$

yields a pointwise nonsingular matrix function
$Q=[\>\Phi\>\>\Phi^{\prime}\>]$. Since $E\dot{\Phi}=A\Phi$, we obtain

$(E,A)\sim([\>E\Phi\>\>E\Phi^{\prime}\>],[\>A\Phi\>\>A\Phi^{\prime}\>]-[\>E\dot{\Phi}\>\>E\dot{\Phi}^{\prime}\>])=([\>E_{1}\>\>E_{2}\>],[\>0\>\>A_{2}\>]).$

In this relation, $E_{1}$ has full column rank $d$. To see this, suppose that
$\mathop{\rm rank}\nolimits E_{1}(\hat{t})<d$ for some
$\hat{t}\in{\mathbb{I}}$. Then there would exist a vector $w\neq 0$ with

$E_{1}(\hat{t})w=0.$

Defining in this situation

$f(t)=\left\\{\begin{array}[]{ll}{1\over t-\hat{t}}E_{1}(t)w&\hbox{for
$t\neq\hat{t}$},\\\\[2.84526pt] {d\over dt}(E_{1}(t)w)&\hbox{for
$t=\hat{t}$},\end{array}\right.$

we would obtain a smooth inhomogeneity $f$. The function $x$ given by

$x(t)=\left[\begin{array}[]{c}\log(|t-\hat{t}|)w\\\ 0\end{array}\right]$

would then solve

$[\>E_{1}(t)\>\>E_{2}(t)\>]\dot{x}=[\>0\>\>A_{2}(t)\>]x+f(t)$

on ${\mathbb{I}}\setminus\\{\hat{t}\\}$ in contradiction to the assumption of
unique solvability, which includes by definition that solutions are defined on
the entire interval ${\mathbb{I}}$.

Hence, since $E_{1}$ has full column rank, there exists a smooth, pointwise
nonsingular matrix function $P$ with

$PE_{1}=\left[\begin{array}[]{c}I_{d}\\\ 0\end{array}\right],$

and thus

$(E,A)\sim\left(\left[\begin{array}[]{cc}I_{d}&E_{12}\\\
0&E_{22}\end{array}\right],\left[\begin{array}[]{cc}0&A_{12}\\\
0&A_{22}\end{array}\right]\right).$

The equation

$E_{22}(t)\dot{x}_{2}=A_{22}(t)x_{2}$

only admits the trivial solution. To see this, suppose that $x_{2}\neq 0$ is a
nontrivial solution and $x_{1}$ a solution of the ODE

$\dot{x}_{1}+E_{12}(t)\dot{x}_{2}(t)=A_{22}(t)x_{2}(t).$

Then we obtain

$[\>E_{1}(t)\>\>E_{2}(t)\>]\left[\begin{array}[]{c}\dot{x}_{1}(t)\\\
\dot{x}_{2}(t)\end{array}\right]=[\>0\>\>A_{2}(t)\>]\left[\begin{array}[]{c}x_{1}(t)\\\
x_{2}(t)\end{array}\right].$

Since $x_{2}\neq 0$, transforming back gives

$E(t)Q(t)\left[\begin{array}[]{c}\dot{x}_{1}(t)\\\
\dot{x}_{2}(t)\end{array}\right]=(A(t)Q(t)-E(t)\dot{Q}(t))\left[\begin{array}[]{c}x_{1}(t)\\\
x_{2}(t)\end{array}\right]$

or $E(t)\dot{x}(t)=A(t)x(t)$ with

$x=Q\left[\begin{array}[]{c}x_{1}\\\ x_{2}\end{array}\right]\neq 0,\quad
x\not\in\mathop{\rm span}\nolimits\\{\phi_{1},\ldots,\phi_{d}\\}.$

But this contradicts the construction of $\phi_{1},\ldots,\phi_{d}$, and hence
(15) holds.

In the presence of symmetries, we of course want to maintain these properties,
which requires to restrict the allowed equivalence transformations. We will
make use of the following notions and properties.

###### Definition 2.6

Two pairs $(E_{i},A_{i})$, $E_{i},A_{i}\in
C({\mathbb{I}},{\mathbb{R}}^{n,n})$, $i=1,2$, of matrix functions are called
congruent if there exist a pointwise nonsingular matrix function $Q\in
C^{1}({\mathbb{I}},{\mathbb{R}}^{n,n})$ such that

$E_{2}=Q^{T}E_{1}Q,\quad A_{2}=Q^{T}A_{1}Q-Q^{T}E_{1}\dot{Q}$ (17)

as equality of functions. We write $(E_{1},A_{1})\equiv(E_{2},A_{2})$.

Again, it is easy to see that the relation defined in Definition 2.6 indeed is
an equivalence relation, see [22].

The following result then modifies Theorem 2.3 provided some symmetry property
holds.

###### Theorem 2.7

Let $E\in C^{k}({\mathbb{I}},{\mathbb{R}}^{n,n})$,
$k\in{\mathbb{N}}_{0}\cup\\{\infty\\}$, with $\mathop{\rm rank}\nolimits
E(t)=r$ for all $t\in{\mathbb{I}}$ and let $\mathop{\rm kernel}\nolimits
E(t)^{T}=\mathop{\rm kernel}\nolimits E(t)$ for all $t\in{\mathbb{I}}$. Then
there exist a pointwise orthogonal function $Q\in
C^{k}({\mathbb{I}},{\mathbb{R}}^{n,n})$ such that

$Q^{T}EQ=\begin{bmatrix}\Sigma&0\\\ 0&0\end{bmatrix}$ (18)

with pointwise nonsingular $\Sigma\in C^{k}({\mathbb{I}},{\mathbb{R}}^{r,r})$.

Proof. With $Q=V$ from Theorem 2.3, the pointwise property $\mathop{\rm
kernel}\nolimits E^{T}=\mathop{\rm kernel}\nolimits E$ allows to choose
$U=Q^{T}$.

In the next two sections we will employ these preliminary results to derive
canonical forms for self-adjoint and skew-adjoint DAEs.

##  3 Canonical forms for self-adjoint pairs of matrix functions

In this section we study canonical forms under congruence for self-adjoint
DAEs. We first show that congruence preserves the self-adjoint structure.

###### Lemma 3.1

Consider two pairs of matrix functions $(E,A)$ and
$(\widetilde{E},\widetilde{A})$ that are congruent and let $(E,A)$ be self-
adjoint. Then $(\widetilde{E},\widetilde{A})$ is self-adjoint as well.

Proof. Let

$\widetilde{E}=Q^{T}EQ,\quad\widetilde{A}=Q^{T}AQ-Q^{T}E\dot{Q}$

and assume that (2) holds. Then

$\widetilde{E}^{T}=Q^{T}E^{T}Q=-Q^{T}EQ=-\widetilde{E}$

and

$\begin{array}[]{l}\widetilde{A}^{T}=Q^{T}A^{T}Q-\dot{Q}^{T}E^{T}Q=Q^{T}(A+\dot{E})Q+\dot{Q}^{T}EQ\\\\[2.84526pt]
\hphantom{\widetilde{A}^{T}}=Q^{T}AQ-Q^{T}E\dot{Q}+Q^{T}E\dot{Q}+Q^{T}\dot{E}Q+\dot{Q}^{T}EQ\\\
\hphantom{\widetilde{A}^{T}}=(Q^{T}AQ-Q^{T}E\dot{Q})+\vphantom{\tilde{E}}{\textstyle{d\over
dt}}(Q^{T}EQ)=\widetilde{A}+\dot{\widetilde{E}}.\end{array}$

For self-adjoint pairs the following local canonical form under pointwise
orthogonal congruence transformations is due to [22] stated here for the
special case of a regular pair of matrix functions.

###### Theorem 3.2

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth and let $(E,A)$ be skew-adjoint. Then there exist pairwise
disjoint open intervals ${\mathbb{I}}_{j}$, $j\in{\mathbb{N}}$, with (11) such
that on every ${\mathbb{I}}_{j}$ there exists a pointwise orthogonal $Q\in
C({\mathbb{I}},{\mathbb{R}}^{n,n})$ with

$(Q^{T}EQ,Q^{T}\\!AQ-Q^{T}E\dot{Q})=\left(\left[\begin{array}[]{c|cc|c}*&*&*&E_{14}\\\
\hline\cr*&\Delta&0&0\\\ &0&0&0\\\ \hline\cr
E_{41}&0&0&0\end{array}\right],\left[\begin{array}[]{c|cc|c}*&*&*&A_{14}\\\
\hline\cr*&\Sigma_{11}&\Sigma_{12}&0\\\ &\Sigma_{21}&\Sigma_{22}&0\\\
\hline\cr A_{41}&0&0&0\end{array}\right]\right),$ (19)

where

$E_{14}=\left[\begin{array}[]{cccc}*&\cdots&*&0\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ 0\end{array}\right],\quad
A_{14}=\left[\begin{array}[]{cccc}*&\cdots&*&\Gamma_{w}\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ \Gamma_{1}\end{array}\right]$ (20)

and $\Delta,\Sigma_{22},\Gamma_{1},\ldots,\Gamma_{w}$ are pointwise
nonsingular. Furthermore,

$\Delta^{T}=-\Delta,\quad\Sigma_{11}^{T}=\Sigma_{11}^{\vphantom{(}}+\dot{\Delta},\quad\Sigma_{21}^{T}=\Sigma_{12}^{\vphantom{(}},\quad\Sigma_{22}^{T}=\Sigma_{22}^{\vphantom{(}},\quad
A_{41}^{T}=A_{14}^{\vphantom{(}}+\dot{E}_{14}.$ (21)

Theorem 3.2 can be further refined by allowing for a restricted class of non-
orthogonal transformations, see again [22].

###### Theorem 3.3

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth and let $(E,A)$ be skew-adjoint. Then there exist pairwise
disjoint open intervals ${\mathbb{I}}_{j}$, $j\in{\mathbb{N}}$, with (11) such
that on every ${\mathbb{I}}_{j}$ there exists a pointwise nonsingular $Q\in
C({\mathbb{I}},{\mathbb{R}}^{n,n})$ with

$(Q^{T}EQ,Q^{T}\\!AQ-Q^{T}E\dot{Q})=\left(\left[\begin{array}[]{c|cc|c}*&*&*&E_{14}\\\
\hline\cr*&J&0&0\\\ &0&0&0\\\ \hline\cr
E_{41}&0&0&0\end{array}\right],\left[\begin{array}[]{c|cc|c}*&*&*&A_{14}\\\
\hline\cr*&C&0&0\\\ &0&\Sigma_{22}&0\\\ \hline\cr
A_{41}&0&0&0\end{array}\right]\right),$ (22)

where

$E_{14}=\left[\begin{array}[]{ccccc}*&\cdots&*&0\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ 0\end{array}\right],\quad
A_{14}=\left[\begin{array}[]{ccccc}*&\cdots&*&I\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ I\end{array}\right].$ (23)

and $\Sigma_{22}$ pointwise nonsingular. Furthermore,

$J=\left[\begin{array}[]{cc}0&I_{p}\\\ -I_{p}&0\end{array}\right],\quad
C^{T}=C^{\vphantom{(}},\quad\Sigma_{22}^{T}=\Sigma_{22}^{\vphantom{(}},\quad
A_{41}^{T}=A_{14}^{\vphantom{(}}+\dot{E}_{14}^{\vphantom{(}}.$ (24)

By successively resolving the algebraic equations in the fourth, third and
first row of the DAE with solution
$[\>x_{1}^{T}\>\>x_{2}^{T}\>\>x_{3}^{T}\>\>x_{4}^{T}\>]^{T}$ and inhomogeneity
$[\>f_{1}^{T}\>\>f_{2}^{T}\>\>f_{3}^{T}\>\>f_{4}^{T}\>]^{T}$ associated with
(22), we can directly solve for $x_{1}$ in terms of linear combinations of
derivatives of $f_{4}$, for $x_{3}$ in terms of $x_{1}$ and $f_{3}$, and for
$x_{4}$ in terms of linear combinations of derivatives of $f_{1}$ and all
other components. The only dynamic behavior related to (22) is described by
the second block row. Inserting $x_{1}$ obtained from the last block row and
calling the updated inhomogeneity $\widetilde{f}_{2}$, the associated ODE
reads

$\dot{x}_{2}=J^{-1}C(t)x_{2}+J^{-1}\widetilde{f}_{2}(t).$ (25)

The matrix function $M=J^{-1}C$ satisfies $M^{T}J-JM=0$ and lies therefore
pointwise in the _Lie algebra of Hamiltonian matrices_ , see e.g. [15]. Thus,
the flow corresponding to (25) defined by

$\dot{\Phi}_{2}=J^{-1}C(t)\Phi_{2},\quad\Phi_{2}(t_{0})=I$ (26)

satisfies $\Phi_{2}^{T}J\Phi_{2}=J$, see e.g. [14], and is therefore
symplectic.

In [22], also a global canonical form for self-adjoint DAEs and the associated
pairs $(E,A)$ was derived. The following modified result differs slightly in
the assumptions, and, moreover, the resulting canonical form is more refined.

###### Theorem 3.4

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth and let $(E,A)$ be self-adjoint. Then we have

$(E,A)\equiv\left(\left[\begin{array}[]{ccc}0&I_{p}&0\\\ -I_{p}&0&0\\\
0&0&E_{33}\end{array}\right],\left[\begin{array}[]{ccc}0&0&0\\\
0&A_{22}&A_{23}\\\ 0&A_{32}&A_{33}\end{array}\right]\right),$ (27)

where

$E_{33}(t)\dot{x}_{3}=A_{33}(t)x_{3}+f_{3}(t),$ (28)

is uniquely solvable for every sufficiently smooth $f_{3}$ without specifying
initial conditions. Furthermore,

$E_{33}^{T}=-E_{33}^{\vphantom{(}},\quad
A_{22}^{T}=A_{22}^{\vphantom{(}},\quad A_{32}^{T}=A_{23}^{\vphantom{(}},\quad
A_{33}^{T}=A_{33}+\dot{E}_{33}.$ (29)

Proof. According to Theorem 3.3, the size $d$ of the differential part is
given by $d=2p$. This implies that the solution space of the homogeneous
equation

$E(t)\dot{x}=A(t)x$

is of dimension $2p$. If $p=0$, then the first two blocks are missing and the
claim holds trivially by assumption.

Let therefore $p\neq 0$ and let $\\{\phi_{1},\ldots,\phi_{2p}\\}$ be a basis
of the solution space. Setting $\Phi=[\>\phi_{1}\>\>\cdots\>\>\phi_{2p}\>]$,
we have

$\mathop{\rm rank}\nolimits\Phi(t)=2p$ for all $t\in{\mathbb{I}}$

as in the general case of Theorem 2.5.

Hence, by Theorem 2.3 there exists a smooth, pointwise nonsingular matrix
function $U$ with

$U^{H}\Phi=\left[\begin{array}[]{c}I_{2p}\\\ 0\end{array}\right].$

Defining

$\Phi^{\prime}=U\left[\begin{array}[]{c}0\\\ I_{a}\end{array}\right]$

with $a=n-2p$ yields a pointwise nonsingular matrix function
$Q=[\>\Phi\>\>\Phi^{\prime}\>]$. Since $E\dot{\Phi}=A\Phi$, we obtain

$(\widetilde{E},\widetilde{A})=(Q^{T}EQ,Q^{T}AQ-Q^{T}E\dot{Q})$

with

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{cc}\Phi^{T}E\Phi&\Phi^{T}E\Phi^{\prime}\\\
\Phi^{\prime T}E\Phi&\Phi^{\prime
T}E\Phi^{\prime}\end{array}\right]=\left[\begin{array}[]{cc}\widetilde{E}_{11}&\widetilde{E}_{12}\\\
-\widetilde{E}_{12}^{T}&\widetilde{E}_{22}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{cc}\Phi^{T}(A\Phi-E\dot{\Phi})&\Phi^{T}(A\Phi^{\prime}-E\dot{\Phi}^{\prime})\\\
\Phi^{\prime T}(A\Phi-E\dot{\Phi})&\Phi^{\prime
T}(A\Phi^{\prime}-E\dot{\Phi}^{\prime})\end{array}\right]=\left[\begin{array}[]{cc}0&\widetilde{A}_{12}\\\
0&\widetilde{A}_{22}\end{array}\right].\end{array}$

To simplify the notation, we omit now and later at similar instances the
tildes thus re-using the same notation for possibly different quantities and
write

$(E,A)\equiv\left(\left[\begin{array}[]{cc}E_{11}&E_{12}\\\
-E_{12}^{T}&E_{22}\end{array}\right],\left[\begin{array}[]{cc}0&A_{12}\\\
0&A_{22}\end{array}\right]\right).$

As in the general case, we can conclude that

$\mathop{\rm rank}\nolimits\left[\begin{array}[]{c}E_{11}\\\
-E_{12}^{T}\end{array}\right]=2p.$

Since self-adjointness is conserved, we additionally have

$E_{11}^{T}=-E_{11},\quad E_{22}^{T}=-E_{22},\quad 0=\dot{E}_{11},\quad
0=A_{12}+\dot{E}_{12},\quad A_{22}^{T}=A_{22}+\dot{E}_{22}.$

In particular, $E_{11}$ is constant and skew-symmetric. Hence, see [22], there
exists an orthogonal symplectic matrix $U\in{\mathbb{R}}^{2p,2p}$ with

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{cc}U^{T}&0\\\
0&I_{a}\end{array}\right]\left[\begin{array}[]{cc}E_{11}&E_{12}\\\
-E_{12}^{T}&E_{22}\end{array}\right]\left[\begin{array}[]{cc}U&0\\\
0&I_{a}\end{array}\right]=\left[\begin{array}[]{ccc}0&\widetilde{E}_{12}&\widetilde{E}_{13}\\\
-\widetilde{E}_{12}^{T}&\widetilde{E}_{22}&\widetilde{E}_{23}\\\
-\widetilde{E}_{13}^{T}&-\widetilde{E}_{23}^{T}&\widetilde{E}_{33}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{cc}U^{T}&0\\\
0&I_{a}\end{array}\right]\left[\begin{array}[]{cc}0&A_{12}\\\
0&A_{22}\end{array}\right]\left[\begin{array}[]{cc}U&0\\\
0&I_{a}\end{array}\right]=\left[\begin{array}[]{ccc}0&0&\widetilde{A}_{13}\\\
0&0&\widetilde{A}_{23}\\\
0&0&\widetilde{A}_{33}\end{array}\right],\end{array}$

where $\widetilde{E}_{12}\in{\mathbb{R}}^{p,p}$.

Omitting again the tildes, we write

$(E,A)\equiv\left(\left[\begin{array}[]{ccc}0&E_{12}&E_{13}\\\
-E_{12}^{T}&E_{22}&E_{23}\\\
-E_{13}^{T}&-E_{23}^{T}&E_{33}\end{array}\right],\left[\begin{array}[]{ccc}0&0&A_{13}\\\
0&0&A_{23}\\\ 0&0&A_{33}\end{array}\right]\right).$

Conservation of self-adjointness and full rank of the leading block yields

$E_{22}^{T}=-E_{22},\quad E_{33}^{T}=-E_{33},\quad 0=\dot{E}_{12},\quad
0=\dot{E}_{22},\quad\mathop{\rm rank}\nolimits[\>E_{12}\>\>E_{13}\>]=p.$

Hence, there exists a smooth, pointwise nonsingular matrix function $V$ with

$[\>E_{12}\>\>E_{13}\>]V=[\>I_{p}\>\>0\>]$

leading to

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{cc}I_{p}&0\\\
0&V^{T}\end{array}\right]\left[\begin{array}[]{ccc}0&E_{12}&E_{13}\\\
-E_{12}^{T}&E_{22}&E_{23}\\\
-E_{13}^{T}&-E_{23}^{T}&E_{33}\end{array}\right]\left[\begin{array}[]{cc}I_{p}&0\\\
0&V\end{array}\right]=\left[\begin{array}[]{ccc}0&I_{p}&0\\\
-I_{p}&\widetilde{E}_{22}&\widetilde{E}_{23}\\\
0&-\widetilde{E}_{23}^{T}&\widetilde{E}_{33}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{cc}I_{p}&0\\\
0&V^{T}\end{array}\right]\left[\begin{array}[]{ccc}0&0&\widetilde{A}_{13}\\\
0&0&\widetilde{A}_{23}\\\
0&0&\widetilde{A}_{33}\end{array}\right]\left[\begin{array}[]{cc}I_{p}&0\\\
0&V\end{array}\right]\\\ \qquad\qquad{}-\left[\begin{array}[]{cc}I_{p}&0\\\
0&V^{T}\end{array}\right]\left[\begin{array}[]{ccc}0&E_{12}&E_{13}\\\
-E_{12}^{T}&E_{22}&E_{23}\\\
-E_{13}^{T}&-E_{23}^{T}&E_{33}\end{array}\right]\left[\begin{array}[]{cc}0&0\\\
0&\dot{V}\end{array}\right]=\left[\begin{array}[]{ccc}0&\widetilde{A}_{12}&\widetilde{A}_{13}\\\
0&\widetilde{A}_{22}&\widetilde{A}_{23}\\\
0&\widetilde{A}_{32}&\widetilde{A}_{33}\end{array}\right].\end{array}$

Omitting the tildes again, we write

$(E,A)\equiv\left(\left[\begin{array}[]{ccc}0&I_{p}&0\\\
-I_{p}&E_{22}&E_{23}\\\
0&-E_{23}^{T}&E_{33}\end{array}\right],\left[\begin{array}[]{ccc}0&A_{12}&A_{13}\\\
0&A_{22}&A_{23}\\\ 0&A_{32}&A_{33}\end{array}\right]\right).$

Conservation of self-adjointness yields

$E_{22}^{T}=-E_{22},\quad E_{33}^{T}=-E_{33},\quad 0=A_{12},\quad 0=A_{13}.$

Finally, after a congruence transformation with,

$\left[\begin{array}[]{ccc}I_{p}&\frac{1}{2}E_{22}&E_{23}\\\ 0&I_{p}&0\\\
0&0&I_{a}\end{array}\right],$

we arrive at

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{ccc}I_{p}&0&0\\\
-\frac{1}{2}E_{22}&I_{p}&0\\\
E_{23}^{T}&0&I_{a}\end{array}\right]\left[\begin{array}[]{ccc}0&I_{p}&0\\\
-I_{p}&E_{22}&E_{23}\\\
0&-E_{23}^{T}&E_{33}\end{array}\right]\left[\begin{array}[]{ccc}I_{p}&\frac{1}{2}E_{22}&E_{23}\\\
0&I_{p}&0\\\
0&0&I_{a}\end{array}\right]=\left[\begin{array}[]{ccc}0&I_{p}&0\\\
-I_{p}&0&0\\\ 0&0&\widetilde{E}_{33}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{ccc}I_{p}&0&0\\\
-\frac{1}{2}E_{22}&I_{p}&0\\\
E_{23}^{T}&0&I_{a}\end{array}\right]\left[\begin{array}[]{ccc}0&A_{12}&A_{13}\\\
0&A_{22}&A_{23}\\\
0&A_{32}&A_{33}\end{array}\right]\left[\begin{array}[]{ccc}I_{p}&\frac{1}{2}E_{22}&E_{23}\\\
0&I_{p}&0\\\ 0&0&I_{a}\end{array}\right]\\\
\qquad\qquad{}-\left[\begin{array}[]{ccc}I_{p}&0&0\\\
-\frac{1}{2}E_{22}&I_{p}&0\\\
E_{23}^{T}&0&I_{a}\end{array}\right]\left[\begin{array}[]{ccc}0&I_{p}&0\\\
-I_{p}&E_{22}&E_{23}\\\
0&-E_{23}^{T}&E_{33}\end{array}\right]\left[\begin{array}[]{ccc}0&\frac{1}{2}\dot{E}_{22}&\dot{E}_{23}\\\
0&0&0\\\ 0&0&0\end{array}\right]=\left[\begin{array}[]{ccc}0&0&0\\\
0&\widetilde{A}_{22}&\widetilde{A}_{23}\\\
0&\widetilde{A}_{32}&\widetilde{A}_{33}\end{array}\right].\end{array}$

Omitting again the tildes then yields

$(E,A)\equiv\left(\left[\begin{array}[]{ccc}0&I_{p}&0\\\ -I_{p}&0&0\\\
0&0&E_{33}\end{array}\right],\left[\begin{array}[]{ccc}0&0&0\\\
0&A_{22}&A_{23}\\\ 0&A_{32}&A_{33}\end{array}\right]\right)$

which is just (27), where (28) follows along the same lines as in the general
case and (29) follows by the conservation of self-adjointness.

After having derived canonical forms for self-adjoint DAEs, in the next
section we derive an analogous form for skew-adjoint DAEs.

##  4 Canonical forms for skew-adjoint pairs of matrix functions

In this section we show that canonical forms under congruence can be derived
also for skew-adjoint DAEs.

###### Lemma 4.1

Consider two pairs of matrix functions $(E,A)$ and
$(\widetilde{E},\widetilde{A})$ that are congruent and let $(E,A)$ be skew-
adjoint. Then $(\widetilde{E},\widetilde{A})$ is skew-adjoint as well.

Proof. Let

$\widetilde{E}=Q^{T}EQ,\quad\widetilde{A}=Q^{T}AQ-Q^{T}E\dot{Q}$

and assume that (3) holds. Then

$\widetilde{E}^{T}=Q^{T}E^{T}Q=Q^{T}EQ=\widetilde{E}$

and

$\begin{array}[]{l}\widetilde{A}^{T}=Q^{T}A^{T}Q-\dot{Q}^{T}E^{T}Q=Q^{T}(-A-\dot{E})Q-\dot{Q}^{T}EQ\\\\[2.84526pt]
\hphantom{\widetilde{A}^{T}}=-Q^{T}AQ+Q^{T}E\dot{Q}-Q^{T}E\dot{Q}-Q^{T}\dot{E}Q-\dot{Q}^{T}EQ\\\
\hphantom{\widetilde{A}^{T}}=-(Q^{T}AQ-Q^{T}E\dot{Q})-\vphantom{\tilde{E}}{\textstyle{d\over
dt}}(Q^{T}EQ)=-\widetilde{A}-\dot{\widetilde{E}}.\end{array}$

The following result on a local canonical form under pointwise orthogonal
congruence transformations is due to [37].

###### Theorem 4.2

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth and let $(E,A)$ be skew-adjoint. Then there exist pairwise
disjoint open intervals ${\mathbb{I}}_{j}$, $j\in{\mathbb{N}}$, with (11) such
that on every ${\mathbb{I}}_{j}$ there exists a pointwise orthogonal $Q\in
C({\mathbb{I}},{\mathbb{R}}^{n,n})$ with

$(Q^{T}EQ,Q^{T}\\!AQ-Q^{T}E\dot{Q})=\left(\left[\begin{array}[]{c|cc|c}*&*&*&E_{14}\\\
\hline\cr*&\Delta&0&0\\\ &0&0&0\\\ \hline\cr
E_{41}&0&0&0\end{array}\right],\left[\begin{array}[]{c|cc|c}*&*&*&A_{14}\\\
\hline\cr*&\Sigma_{11}&\Sigma_{12}&0\\\ &\Sigma_{21}&\Sigma_{22}&0\\\
\hline\cr A_{41}&0&0&0\end{array}\right]\right),$ (30)

where

$E_{14}=\left[\begin{array}[]{cccc}*&\cdots&*&0\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ 0\end{array}\right],\quad
A_{14}=\left[\begin{array}[]{cccc}*&\cdots&*&\Gamma_{w}\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ \Gamma_{1}\end{array}\right]$ (31)

and $\Delta,\Sigma_{22},\Gamma_{1},\ldots,\Gamma_{w}$ are pointwise
nonsingular. Furthermore,

$\Delta^{T}=\Delta,\quad\Sigma_{11}^{T}=-\Sigma_{11}^{\vphantom{(}}-\dot{\Delta},\quad\Sigma_{21}^{T}=-\Sigma_{12}^{\vphantom{(}},\quad\Sigma_{22}^{T}=-\Sigma_{22}^{\vphantom{(}},\quad
A_{41}^{T}=-A_{14}^{\vphantom{(}}-\dot{E}_{14}.$ (32)

Theorem 3.2 can be further refined by allowing for a restricted class of non-
orthogonal transformations, which yields the following local canonical form
which, using a recent result of [19], is more refined than that in [37].

###### Theorem 4.3

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth and let $(E,A)$ be skew-adjoint. Then there exist pairwise
disjoint open intervals ${\mathbb{I}}_{j}$, $j\in{\mathbb{N}}$, with (11) such
that on every ${\mathbb{I}}_{j}$ there exists a pointwise nonsingular $Q\in
C({\mathbb{I}},{\mathbb{R}}^{n,n})$ with

$(Q^{T}EQ,Q^{T}\\!AQ-Q^{T}E\dot{Q})=\left(\left[\begin{array}[]{c|cc|c}*&*&*&E_{14}\\\
\hline\cr*&S&0&0\\\ &0&0&0\\\ \hline\cr
E_{41}&0&0&0\end{array}\right],\left[\begin{array}[]{c|cc|c}*&*&*&A_{14}\\\
\hline\cr*&J&0&0\\\ &0&\Sigma_{22}&0\\\ \hline\cr
A_{41}&0&0&0\end{array}\right]\right),$ (33)

where

$E_{14}=\left[\begin{array}[]{ccccc}*&\cdots&*&0\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ 0\end{array}\right],\quad
A_{14}=\left[\begin{array}[]{ccccc}*&\cdots&*&I\\\ \vdots&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}&\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ &\mbox{
\begin{picture}(8.0,7.0)\put(-4.0,1.0){.}\put(0.0,4.0){.}
\put(4.0,7.0){.}\end{picture}}\\\ I\end{array}\right].$ (34)

and $\Sigma_{22}$ pointwise nonsingular. Furthermore,

$S=\left[\begin{array}[]{cc}I_{p}&0\\\ 0&-I_{q}\end{array}\right],\quad
J^{T}=-J^{\vphantom{(}},\quad\Sigma_{22}^{T}=-\Sigma_{22}^{\vphantom{(}},\quad
A_{41}^{T}=-A_{14}^{\vphantom{(}}-\dot{E}_{14}^{\vphantom{(}}.$ (35)

Proof. Compared with the result in [37], to obtain the stated local canonical
form we need a smooth version of Sylvester’s law of inertia, which is proved
in [19]. In particular, this result shows the existence of a smooth
transformation $W$ with $W^{T}\\!\Delta\,W=S$, where $\Delta$ is pointwise
nonsingular and symmetric.

By successively resolving the algebraic equations in the fourth, third and
first row of the DAE with solution
$[\>x_{1}^{T}\>\>x_{2}^{T}\>\>x_{3}^{T}\>\>x_{4}^{T}\>]^{T}$ and inhomogeneity
$[\>f_{1}^{T}\>\>f_{2}^{T}\>\>f_{3}^{T}\>\>f_{4}^{T}\>]^{T}$ associated with
(22), we can directly solve for $x_{1}$ in terms of linear combinations of
derivatives of $f_{4}$, for $x_{3}$ in terms of $x_{1}$ and $f_{3}$, and for
$x_{4}$ in terms of linear combinations of derivatives of $f_{1}$ and all
other components. The only dynamic behavior related to (22) is described by
the second block row. Inserting $x_{1}$ obtained from the last block row and
calling the updated inhomogeneity $\widetilde{f}_{2}$, the associated ODE
reads

$\dot{x}_{2}=S^{-1}J(t)x_{2}+S^{-1}\widetilde{f}_{2}(t).$ (36)

The matrix function $M=S^{-1}J$ has the property that $SM$ is skew-symmetric,
i.e., that $(SM)^{T}=-(SM)$ or $M^{T}S+SM=0$, and lies therefore pointwise in
the Lie algebra belonging to the quadratic Lie group

$O(p,q)=\\{\Phi\in{\mathbb{R}}^{n,n}\mid\Phi^{T}S\Phi=S\\},$

the so-called _generalized orthogonal group with inertia $(p,q,0)$_, see e.g.
[15], implying that the flow belonging to (36) lies in $O(p,q)$. Observe the
special case $p=0$ or $q=0$ where the flow is then orthogonal.

We also obtain a global canonical form.

###### Theorem 4.4

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth and let $(E,A)$ be skew-adjoint. Then we have

$(E,A)\equiv\left(\left[\begin{array}[]{ccc}I_{p}&0&0\\\ 0&-I_{q}&0\\\
0&0&E_{33}\end{array}\right],\left[\begin{array}[]{ccc}0&0&0\\\ 0&0&0\\\
0&0&A_{33}\end{array}\right]\right),$ (37)

where

$E_{33}(t)\dot{x}_{3}=A_{33}(t)x_{3}+f_{3}(t)$ (38)

is uniquely solvable for every sufficiently smooth $f_{3}$ without specifying
initial conditions. Furthermore,

$E_{33}^{T}=E_{33}^{\vphantom{(}},\quad
A_{33}^{T}=-A_{33}-\dot{E}_{33}^{\vphantom{(}}$ (39)

Proof. According to Theorem 4.3, the size $d$ of the differential part is
given by $d=p+q$. This implies that the solution space of the homogeneous
equation

$E(t)\dot{x}=A(t)x$

is of dimension $p+q$. If $p+q=0$, then the first two blocks are missing and
the claim holds trivially by assumption.

Let therefore $p+q\neq 0$ and let $\\{\phi_{1},\ldots,\phi_{p+q}\\}$ be a
basis of the solution space. Setting
$\Phi=[\>\phi_{1}\>\>\cdots\>\>\phi_{p+q}\>]$, we have

$\mathop{\rm rank}\nolimits\Phi(t)=p+q$ for all $t\in{\mathbb{I}}$

as in the general case of Theorem 2.5.

Hence, by Theorem 2.3 there exists a smooth, pointwise nonsingular matrix
function $U$ with

$U^{H}\Phi=\left[\begin{array}[]{c}I_{p+q}\\\ 0\end{array}\right].$

Defining

$\Phi^{\prime}=U\left[\begin{array}[]{c}0\\\ I_{a}\end{array}\right]$

with $a=n-(p+q)$ yields a pointwise nonsingular matrix function
$Q=[\>\Phi\>\>\Phi^{\prime}\>]$. Since $E\dot{\Phi}=A\Phi$, we obtain

$(\widetilde{E},\widetilde{A})=(Q^{T}EQ,Q^{T}AQ-Q^{T}E\dot{Q})$

with

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{cc}\Phi^{T}E\Phi&\Phi^{T}E\Phi^{\prime}\\\
\Phi^{\prime T}E\Phi&\Phi^{\prime
T}E\Phi^{\prime}\end{array}\right]=\left[\begin{array}[]{cc}\widetilde{E}_{11}&\widetilde{E}_{12}\\\
\widetilde{E}_{12}^{T}&\widetilde{E}_{22}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{cc}\Phi^{T}(A\Phi-E\dot{\Phi})&\Phi^{T}(A\Phi^{\prime}-E\dot{\Phi}^{\prime})\\\
\Phi^{\prime T}(A\Phi-E\dot{\Phi})&\Phi^{\prime
T}(A\Phi^{\prime}-E\dot{\Phi}^{\prime})\end{array}\right]=\left[\begin{array}[]{cc}0&\widetilde{A}_{12}\\\
0&\widetilde{A}_{22}\end{array}\right].\end{array}$

As in the proof of Theorem 3.4, we omit now and later at similar instances the
tildes thus re-using the same notation for possibly different quantities and
write

$(E,A)\equiv\left(\left[\begin{array}[]{cc}E_{11}&E_{12}\\\
E_{12}^{T}&E_{22}\end{array}\right],\left[\begin{array}[]{cc}0&A_{12}\\\
0&A_{22}\end{array}\right]\right).$

As in the general case, we can conclude that

$\mathop{\rm rank}\nolimits\left[\begin{array}[]{c}E_{11}\\\
E_{12}^{T}\end{array}\right]=p+q.$

Since skew-adjointness is conserved, we additionally have

$E_{11}^{T}=E_{11},\quad E_{22}^{T}=E_{22},\quad 0=\dot{E}_{11},\quad
0=-A_{12}-\dot{E}_{12},\quad A_{22}^{T}=-A_{22}-\dot{E}_{22}.$

In particular, $E_{11}$ is constant and symmetric. Moreover, in the following
we will show that $E_{11}$ is nonsingular.

Due to Sylvester’s law of inertia, there is a nonsingular matrix
$U\in{\mathbb{R}}^{d,d}$, $d=p+q+r$, with

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{cc}U^{T}&0\\\
0&I_{a}\end{array}\right]\left[\begin{array}[]{cc}E_{11}&E_{12}\\\
E_{12}^{T}&E_{22}\end{array}\right]\left[\begin{array}[]{cc}U&0\\\
0&I_{a}\end{array}\right]=\left[\begin{array}[]{ccc}S&0&\widetilde{E}_{13}\\\
0&0&\widetilde{E}_{23}\\\
\widetilde{E}_{13}^{T}&\widetilde{E}_{23}^{T}&\widetilde{E}_{33}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{cc}U^{T}&0\\\
0&I_{a}\end{array}\right]\left[\begin{array}[]{cc}0&A_{12}\\\
0&A_{22}\end{array}\right]\left[\begin{array}[]{cc}U&0\\\
0&I_{a}\end{array}\right]=\left[\begin{array}[]{ccc}0&0&\widetilde{A}_{13}\\\
0&0&\widetilde{A}_{23}\\\
0&0&\widetilde{A}_{33}\end{array}\right].\end{array}$

where $S=\mathop{\rm diag}\nolimits(I_{p},-I_{q})$.

Omitting again the tildes, we write

$(E,A)\equiv\left(\left[\begin{array}[]{ccc}S&0&E_{13}\\\ 0&0&E_{23}\\\
E_{13}^{T}&E_{23}^{T}&E_{33}\end{array}\right],\left[\begin{array}[]{ccc}0&0&A_{13}\\\
0&0&A_{23}\\\ 0&0&A_{33}\end{array}\right]\right).$

Conservation of skew-adjointness and full rank of the leading block yields

$E_{33}=E_{33}^{T},\quad 0=-A_{13}-\dot{E}_{13},\quad
0=-A_{23}-\dot{E}_{23},\quad A_{33}^{T}=-A_{33}-\dot{E}_{33},\quad\mathop{\rm
rank}\nolimits E_{23}=r.$

Hence, there exists a smooth, pointwise nonsingular matrix function $V$ with

$E_{23}V=[\>I_{r}\>\>0\>]$

leading to

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{ccc}I_{p+q}&0&0\\\
0&I_{r}&0\\\
0&0&V^{T}\end{array}\right]\left[\begin{array}[]{ccc}S&0&E_{13}\\\
0&0&E_{23}\\\
E_{13}^{T}&E_{23}^{T}&E_{33}\end{array}\right]\left[\begin{array}[]{ccc}I_{p+q}&0&0\\\
0&I_{r}&0\\\
0&0&V\end{array}\right]=\left[\begin{array}[]{cccc}S&0&\widetilde{E}_{13}&\widetilde{E}_{14}\\\
0&0&I_{r}&0\\\
\widetilde{E}_{13}^{T}&I_{r}&\widetilde{E}_{33}&\widetilde{E}_{34}\\\
\widetilde{E}_{14}^{T}&0&\widetilde{E}_{34}^{T}&\widetilde{E}_{44}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{ccc}I_{p+q}&0&0\\\ 0&I_{r}&0\\\
0&0&V^{T}\end{array}\right]\left[\begin{array}[]{ccc}0&0&A_{13}\\\
0&0&A_{23}\\\
0&0&A_{33}\end{array}\right]\left[\begin{array}[]{ccc}I_{p+q}&0&0\\\
0&I_{r}&0\\\ 0&0&V\end{array}\right]\\\
\qquad\qquad{}-\left[\begin{array}[]{ccc}I_{p+q}&0&0\\\ 0&I_{r}&0\\\
0&0&V^{T}\end{array}\right]\left[\begin{array}[]{ccc}S&0&E_{13}\\\
0&0&E_{23}\\\
E_{13}^{T}&E_{23}^{T}&E_{33}\end{array}\right]\left[\begin{array}[]{ccc}0&0&0\\\
0&0&0\\\
0&0&\dot{V}\end{array}\right]=\left[\begin{array}[]{cccc}0&0&\widetilde{A}_{13}&\widetilde{A}_{14}\\\
0&0&\widetilde{A}_{23}&\widetilde{A}_{24}\\\
0&0&\widetilde{A}_{33}&\widetilde{A}_{34}\\\
0&0&\widetilde{A}_{43}&\widetilde{A}_{44}\end{array}\right].\end{array}$

Omitting again the tildes, we write

$(E,A)\equiv\left(\left[\begin{array}[]{cccc}S&0&E_{13}&E_{14}\\\
0&0&I_{r}&0\\\ E_{13}^{T}&I_{r}&E_{33}&E_{34}\\\
E_{14}^{T}&0&E_{34}^{T}&E_{44}\end{array}\right],\left[\begin{array}[]{cccc}0&0&A_{13}&A_{14}\\\
0&0&A_{23}&A_{24}\\\ 0&0&A_{33}&A_{34}\\\
0&0&A_{43}&A_{44}\end{array}\right].\right).$

Conservation of self-adjointness yields

$E_{33}^{T}=E_{33},\quad E_{44}^{T}=E_{44},\quad 0=-A_{13}-\dot{E}_{13},\quad
0=-A_{14}-\dot{E}_{14},\quad 0=A_{23},\quad 0=A_{24},$

and

$A_{33}^{T}=-A_{33}-\dot{E}_{33},\quad A_{43}^{T}=-A_{34}-\dot{E}_{34},\quad
A_{44}^{T}=-A_{44}-\dot{E}_{44}.$

Finally, after a congruence transformation with,

$\left[\begin{array}[]{cccc}I_{p+q}&0&-S^{-1}E_{13}&-S^{-1}E_{14}\\\
0&I_{r}&0&0\\\ 0&0&I_{r}&0\\\ 0&0&0&I_{a}\end{array}\right],$

we arrive at

$\begin{array}[]{l}\widetilde{E}=\left[\begin{array}[]{cccc}I_{p+q}&0&0&0\\\
0&I_{r}&0&0\\\ -E_{13}^{T}S^{-1}&0&I_{r}&0\\\
-E_{14}^{T}S^{-1}&0&0&I_{a}\end{array}\right]\left[\begin{array}[]{cccc}S&0&E_{13}&E_{14}\\\
0&0&I_{r}&0\\\ E_{13}^{T}&I_{r}&E_{33}&E_{34}\\\
E_{14}^{T}&0&E_{34}^{T}&E_{44}\end{array}\right]\left[\begin{array}[]{cccc}I_{p+q}&0&-S^{-1}E_{13}&-S^{-1}E_{14}\\\
0&I_{r}&0&0\\\ 0&0&I_{r}&0\\\ 0&0&0&I_{a}\end{array}\right]\\\
\qquad\qquad\qquad\qquad{}=\left[\begin{array}[]{cccc}S&0&0&0\\\
0&0&I_{r}&0\\\ 0&I_{r}&\widetilde{E}_{33}&\widetilde{E}_{34}\\\
0&0&\widetilde{E}_{34}^{T}&\widetilde{E}_{44}\end{array}\right],\\\
\widetilde{A}=\left[\begin{array}[]{cccc}I_{p+q}&0&0&0\\\ 0&I_{r}&0&0\\\
-E_{13}^{T}S^{-1}&0&I_{r}&0\\\
-E_{14}^{T}S^{-1}&0&0&I_{a}\end{array}\right]\left[\begin{array}[]{cccc}0&0&A_{13}&A_{14}\\\
0&0&A_{23}&A_{24}\\\ 0&0&A_{33}&A_{34}\\\
0&0&A_{43}&A_{44}\end{array}\right]\left[\begin{array}[]{cccc}I_{p+q}&0&-S^{-1}E_{13}&-S^{-1}E_{14}\\\
0&I_{r}&0&0\\\ 0&0&I_{r}&0\\\ 0&0&0&I_{a}\end{array}\right]\\\
\qquad\qquad{}-\left[\begin{array}[]{cccc}I_{p+q}&0&0&0\\\ 0&I_{r}&0&0\\\
-E_{13}^{T}S^{-1}&0&I_{r}&0\\\
-E_{14}^{T}S^{-1}&0&0&I_{a}\end{array}\right]\left[\begin{array}[]{cccc}S&0&E_{13}&E_{14}\\\
0&0&I_{r}&0\\\ E_{13}^{T}&I_{r}&E_{33}&E_{34}\\\
E_{14}^{T}&0&E_{34}^{T}&E_{44}\end{array}\right]\left[\begin{array}[]{cccc}0&0&-S^{-1}\dot{E}_{13}&-S^{-1}\dot{E}_{14}\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\end{array}\right]\\\
\qquad\qquad\qquad\qquad{}=\left[\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\
0&0&\widetilde{A}_{33}&\widetilde{A}_{34}\\\
0&0&\widetilde{A}_{43}&\widetilde{A}_{44}\end{array}\right].\end{array}$

The corresponding homogeneous DAE
$\widetilde{E}\dot{\widetilde{x}}=\widetilde{A}{\widetilde{x}}$ has the form

$\begin{array}[]{l}S\dot{\widetilde{x}}_{1}=0,\\\
\dot{\widetilde{x}}_{3}=0,\\\
\dot{\widetilde{x}}_{2}+\widetilde{E}_{33}\dot{\widetilde{x}}_{3}+\widetilde{E}_{34}\dot{\widetilde{x}}_{4}=\widetilde{A}_{33}{\widetilde{x}}_{3}+\widetilde{A}_{34}{\widetilde{x}}_{4},\\\
\widetilde{E}_{34}^{T}\dot{\widetilde{x}}_{3}+\widetilde{E}_{44}\dot{\widetilde{x}}_{4}=\widetilde{A}_{43}{\widetilde{x}}_{3}+\widetilde{A}_{44}{\widetilde{x}}_{4}.\end{array}$

If we take ${\widetilde{x}}_{1}$ and ${\widetilde{x}}_{3}$ as solutions of the
first two equations, the fourth equation must determine ${\widetilde{x}}_{4}$,
possibly imposing a suitable initial condition. Finally, the third equation
then fixes ${\widetilde{x}}_{2}$ using a suitable initial condition. Hence the
dimension of the solution space is at least $p+q+2r$. But the dimension was
assumed to be $p+q$ implying that $r=0$.

Omitting the tildes again, skipping the second and third block row and column,
which are of zero dimension, and renumbering the indices gives

$(E,A)\equiv\left(\left[\begin{array}[]{ccc}I_{p}&0&0\\\ 0&-I_{q}&0\\\
0&0&E_{33}\end{array}\right],\left[\begin{array}[]{ccc}0&0&0\\\ 0&0&0\\\
0&0&A_{33}\end{array}\right]\right)$

which is just (37), where (38) follows along the same lines as in the general
case and (39) follows by the conservation of skew-adjointness.

The presented canonical forms have some direct consequences in the case that
$E$ is pointwise positive semidefinite.

###### Corollary 4.5

Let $(E,A)$ be regular with $E,A\in C({\mathbb{I}},{\mathbb{R}}^{n,n})$
sufficiently smooth and let $(E,A)$ be skew-adjoint with $E$ pointwise
positive semidefinite. Then in the canonical form (33) $E_{41}=0$ and the flow
associated with the dynamical part of the system is orthogonal.

Proof. Since $Q^{T}EQ$ is positive semidefinite for positive semidefinite $E$,
it follows that the blocks in positions $(1,3)$, $(3,1)$, $(4,1)$ and $(1,4)$
of (33) are zero and that $S$ is pointwise positive definite. Thus $q=0$ and
the flow corresponding to (36) is orthogonal.

###### Remark 4.6

Corollary 4.5 immediately implies also an upper bound on the so-called
differentiation index of the DAE, see [20], which is at most two. This follows
directly from (33), since with $E_{41}=0$ at most one differentiation of the
inhomogeneity is needed.

###### Example 4.7

Consider the circuit in Example 1.3. Reordering the rows and columns in the
order second, third, first, fourth, and fifth, we obtain the skew-adjoint
system in canonical form

$\left[\begin{array}[]{cc|c|cc}C_{1}&0&0&0&0\\\ 0&C_{2}&0&0&0\\\ \hline\cr
0&0&L&0&0\\\ \hline\cr 0&0&0&0&0\\\
0&0&0&0&0\end{array}\right]\left[\begin{array}[]{c}\dot{V}_{1}\\\
\dot{V}_{2}\\\ \hline\cr\dot{I}\\\ \hline\cr\dot{I}_{G}\\\
\dot{I}_{R}\end{array}\right]=\left[\begin{array}[]{cc|c|cc}0&1&0&-1&0\\\
-1&0&1&0&-1\\\ \hline\cr 0&-1&-R_{L}&0&0\\\ \hline\cr 1&0&0&-R_{G}&0\\\
0&1&0&0&-R_{R}\end{array}\right]\left[\begin{array}[]{c}V_{1}\\\ V_{2}\\\
\hline\cr I\\\ \hline\cr I_{G}\\\
I_{R}\end{array}\right]+\left[\begin{array}[]{c}0\\\ 0\\\ \hline\cr 0\\\
\hline\cr 1\\\ 0\end{array}\right]u.$

If $u$ is given and the resistive terms are neglected, i.e.
$R_{L},R_{G},R_{R}=0$, then this is skew-adjoint DAE, in which the last two
equations can be solved as $V_{2}=0$ and $V_{1}=-u$ and from the first two
equations we get $I_{G}=-C_{1}\dot{V}_{1}=C_{1}\dot{u}$ and
$I_{R}=-C_{2}\dot{V}_{2}=0$ and the dynamic equation is

$L\dot{I}=-V_{2}=0,$

which has the orthogonal flow $I=1$.

###### Example 4.8

Consider the linear time-varying DAE (7) in Example 1.4 with $A_{H}=0$ and
$C=0$. Performing a full rank decomposition

$U^{T}B=\left[\begin{array}[]{c}B_{1}\\\ 0\end{array}\right],$

with $B_{1}$ nonsingular (which corresponds to the partitioning of the
velocity into the parts that have divergence zero and nonzero), and applying a
congruence transformation with

$\left[\begin{array}[]{cc}U&0\\\ 0&I\end{array}\right]$

yields a transformed system

$\left[\begin{array}[]{ccc}M_{11}&M_{12}&0\\\ M_{21}&M_{22}&0\\\
0&0&0\end{array}\right]\left[\begin{array}[]{c}\dot{v}_{1}\\\ \dot{v}_{2}\\\
\dot{p}\end{array}\right]=\left[\begin{array}[]{ccc}J_{11}(t)&J_{12}(t)&-B_{1}\\\
J_{21}(t)&J_{22}(t)&0\\\
B_{1}^{T}&0&0\end{array}\right]\left[\begin{array}[]{c}v_{1}\\\ v_{2}\\\
p\end{array}\right]+\left[\begin{array}[]{c}f_{1}(t)\\\ f_{2}(t)\\\
0\end{array}\right].$

The third equation yields $v_{1}=0$ and the first equation gives

$p=B_{1}^{-1}(J_{12}(t)v_{2}-M_{12}\dot{v}_{2}+f_{1}(t)),$

while the underlying dynamics of the system is described by the skew-adjoint
DAE

$M_{22}\dot{v}_{2}=J_{22}(t)v_{2}+f_{2}(t)$

with constant $M_{22}=M_{22}^{T}>0$ and pointwise skew-symmetric $J_{22}$.
After a change of basis with the positive definite square root
$\smash{M_{22}^{1/2}}$ of $M_{22}$ according to
$\smash{\tilde{v}_{2}=M_{22}^{1/2}v_{2}}$ and scaling the equation by its
inverse $\smash{M_{22}^{-1/2}}$, we obtain an ODE system

$\dot{\tilde{v}}_{2}=\tilde{J}_{22}(t)\tilde{v}_{2}+\tilde{f}_{2}(t),$

with pointwise skew-symmetric $\tilde{J}_{22}$, which has an orthogonal flow.

##  5 Conclusions

We have derived local and global canonical forms under congruence
transformations for self-adjoint and skew-adjoint systems of linear variable
coefficient differential-algebraic equations. The associated flows for the
dynamical part of the system are shown to be symplectic or in the generalized
orthogonal groups. The results are illustrated at the hand of examples from
electrical network and flow simulation.

## Acknowledgement

Volker Mehrmann was partially supported by Deutsche Forschungsgemeinschaft
through the Excellence Cluster Math+ in Berlin and Priority Program 1984
‘Hybride und multimodale Energiesysteme: Systemtheoretische Methoden für die
Transformation und den Betrieb komplexer Netze’.

## References

  * [1] R. Altmann, V. Mehrmann, and B. Unger.  Port-Hamiltonian formulations of poroelastic network models.  Math. and Comp. Modelling Dynamical Systems, 27:429–452, 2021.  http://arxiv.org/abs/2012.01949.
  * [2] R. Altmann and P. Schulze.  On the port-Hamiltonian structure of the Navier-Stokes equations for reactive flows.  Systems & Control Letters, 100:51–55, 2017.
  * [3] C. Beattie, V. Mehrmann, H. Xu, and H. Zwart.  Port-Hamiltonian descriptor systems.  Math. Control, Signals, Sys., 30(17):1–27, 2018.
  * [4] F. Brezzi and M. Fortin.  Mixed and Hybrid Finite Element Methods.  Springer-Verlag, New York, NY, 1991.
  * [5] S. L. Campbell.  A general form for solvable linear time varying singular systems of differential equations.  SIAM J. Math. Anal., 18:1101–1115, 1987.
  * [6] V. Doležal.  The existence of a continuous basis of a certain subspace of $E_{r}$ which depends on a parameter.  Cas. Pro. Pest. Mat., 89:466–468, 1964.
  * [7] P. Domschke, B. Hiller, J. Lang, V. Mehrmann, R. Morandin, and C. Tischendorf.  Gas network modeling: An overview.  Preprint, Collaborative Research Center TRR 154, 2021.
  * [8] H. Egger.  Structure preserving approximation of dissipative evolution problems.  Numer. Math., 143(1):85–106, 2019.
  * [9] E. Eich-Soellner and C. Führer.  Numerical Methods in Multibody Systems.  Teubner Verlag, Stuttgart, Germany, 1998.
  * [10] E. Emmrich and V. Mehrmann.  Operator differential-algebraic equations arising in fluid dynamics.  Comp. Meth. in Appl. Math., 13(4):443–470, 2013.
  * [11] L Formaggia, J.-F. Gerbeau, F. Nobile, and A. Quarteroni.  On the coupling of 3D and 1D Navier–Stokes equations for flow problems in compliant vessels.  Comp. Meth. Appl. Mech. Eng., 191(6-7):561–582, 2001.
  * [12] L. Formaggia, A. Quarteroni, and A. Veneziani.  Cardiovascular Mathematics: Modeling and simulation of the circulatory system, volume 1.  Springer Science & Business Media, 2010.
  * [13] G. Golo, A. J. van der Schaft, P. C. Breedveld, and B. M. Maschke.  Hamiltonian formulation of bond graphs.  In R. Johansson and A. Rantzer, editors, Nonlinear and Hybrid Systems in Automotive Control, pages 351–372. Springer, Heidelberg, 2003.
  * [14] E. Hairer, C. Lubich, and G. Wanner.  Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations.  Springer-Verlag, Berlin, Germany, 2002.
  * [15] B. C. Hall.  Lie groups, Lie algebras, and representations: an elementary introduction, volume 10.  Springer, 2003.
  * [16] S.-A. Hauschild, N. Marheineke, V. Mehrmann, J. Mohring, A. Moses Badlyan, M. Rein, and M. Schmidt.  Port-Hamiltonian modeling of disctrict heating networks.  DAE Forum, Porogress in Differential-Algebraic Equations II, 2020\.  333–355.
  * [17] M. Hiller and K. Hirsch.  Multibody system dynamics and mechatronics.  Z. Angew. Math. Mech., 86(2):87–109, 2006.
  * [18] B. Jacob and H. Zwart.  Linear port-Hamiltonian systems on infinite-dimensional spaces.  Operator Theory: Advances and Applications, 223. Birkhäuser/Springer Basel AG, Basel CH, 2012.
  * [19] P. Kunkel.  A smooth version of Sylvester’s law of inertia and its numerical realization.  Electr. Trans. Num. Anal., 36:542–560, 2020.
  * [20] P. Kunkel and V. Mehrmann.  Differential-Algebraic Equations. Analysis and Numerical Solution.  EMS Publishing House, Zürich, Switzerland, 2006.
  * [21] P. Kunkel and V. Mehrmann.  Optimal control for unstructured nonlinear differential-algebraic equations of arbitrary index.  Math. Control, Signals, Sys., 20:227–269, 2008.
  * [22] P. Kunkel, V. Mehrmann, and L. Scholz.  Self-adjoint differential-algebraic equations.  Math. Control Signals Syst., 26:47–76, 2014.
  * [23] W. Layton.  Introduction to the Numerical Analysis of Incompressible Viscous Flows.  SIAM Publications, Philadelphia, PA, 2008.
  * [24] B.J. Leimkuhler and S. Reich.  Symplectic integration of constrained Hamiltonian systems.  Mathematics of Computation, 63(208):589–605, 1994.
  * [25] V. Mehrmann and R. Morandin.  Structure-preserving discretization for port-Hamiltonian descriptor systems.  In 58th IEEE Conference on Decision and Control (CDC), Nice, pages 6863–6868. IEEE, 2019.
  * [26] R. Ortega, A. J. van der Schaft, Y. Mareels, and B. M. Maschke.  Putting energy back in control.  Control Syst. Mag., 21:18–33, 2001.
  * [27] A. Quarteroni and L. Formaggia.  Mathematical modelling and numerical simulation of the cardiovascular system.  Handbook of numerical analysis, 12:3–127, 2004.
  * [28] A. Quarteroni, A. Manzoni, and C. Vergara.  The cardiovascular system: mathematical modelling, numerical algorithms and clinical applications.  Acta Numerica, 26:365–590, 2017.
  * [29] A. Quarteroni and S. Quarteroni.  Numerical models for differential problems, volume 2.  Springer, 2009.
  * [30] T. Quarteroni, A.and Lassila, S. Rossi, and R. Ruiz-Baier.  Integrated heart—coupling multiscale and multiphysics models for the simulation of the cardiac function.  Comp. Meth. Appl. Mech. Eng., 314:345–407, 2017.
  * [31] R. Rannacher.  Finite element methods for the incompressible Navier-Stokes equations.  In P. Galdi, J. Heywood, and R. Rannacher, editors, Fundamental Directions in Mathematical Fluid Mechanics, pages 191–293. Birkhäuser-Verlag, Basel, Switzerland, 2000.
  * [32] H.-G. Roos, M. Stynes, and L. Tobiska.  Robust Numerical Methods for Singularly Perturbed Differential Equations.  Springer Science & Business Media, 2008.
  * [33] A. J. van der Schaft.  Port-Hamiltonian systems: network modeling and control of nonlinear physical systems.  In Advanced Dynamics and Control of Structures and Machines, CISM Courses and Lectures, Vol. 444. Springer Verlag, New York, N.Y., 2004.
  * [34] A. J. van der Schaft.  Port-Hamiltonian systems: an introductory survey.  In J. L. Verona M. Sanz-Sole and J. Verdura, editors, Proc. of the International Congress of Mathematicians, vol. III, Invited Lectures, pages 1339–1365, Madrid, Spain, 2006.
  * [35] W. Schiehlen.  Advanced multibody system dynamics.  Kluwer Academic Publishers, Stuttgart, Germany, 1993.
  * [36] K. Schlacher and A. Kugi.  Automatic control of mechatronic systems.  Int. J. Appl. Math. Comput. Sci., 11(1):131–164, 2001.
  * [37] L. Scholz.  Condensed forms for linear port-Hamiltonian descriptor systems.  Electr. J. Lin. Alg., 35:65–89, 2019.
  * [38] R. Temam.  Navier-Stokes equations. Theory and numerical analysis.  North Holland, Amsterdam, The Netherlands, 1977.
  * [39] A. J. van der Schaft.  Port-Hamiltonian differential-algebraic systems.  In Surveys in Differential-Algebraic Equations. I, pages 173–226. Springer-Verlag, Heidelberg, 2013.
  * [40] A. J. van der Schaft and D. Jeltsema.  Port-Hamiltonian systems theory: An introductory overview.  Foundations and Trends in Systems and Control, 1(2-3):173–378, 2014\.


# The role of sign indefinite invariants in shaping turbulent cascades

Michal Shavit, Oliver Bühler and Jalal Shatah  Courant Institute of
Mathematical Sciences, New York University, NY 10012, USA.

###### Abstract

We highlight a non-canonical yet natural choice of variables for an efficient
derivation of a kinetic equation for the energy density in non-isotropic
systems, including internal gravity waves on a vertical plane, inertial and
Rossby waves. The existence of a second quadratic invariant simplifies the
kinetic equation and leads to extra conservation laws for resonant
interactions. We analytically determine the scaling of the radial turbulent
energy spectrum. Our findings suggest the existence of an inverse energy
cascade of internal gravity waves, from small to large scales, in practically
relevant scenarios.

††preprint: APS/123-QED

Introduction Strongly dispersive waves are ubiquitous in geophysical fluid
dynamics, where they occur on scales from centimeters to thousands of
kilometres and contribute in an essential and intricate way to the long-term
nonlinear dynamics of the climate system [1, 2, 3, 4]. Examples include
surface waves, internal inertia–gravity waves, and Rossby waves, all of which
owe their existence to some combination of gravity, rotation, and curvature of
the Earth. Many of these waves are far too small in scale to be resolvable
numerically, making their study a pressing issue for theoretical modeling and
investigation. For small-amplitude waves the methods of wave turbulence theory
can play an important part in this, because they produce a closed kinetic
equation for the slow evolution of the averaged spectral energy density. There
has been significant progress for idealized model systems [e.g. 5, 6], but so
far this has not yet been translated to systems of direct geophysical
interest. Arguably, progress has been hampered by the extremely cumbersome
form taken by the relevant equations when attempting to shoe-horn them into
classical wave turbulence theory, which was formulated in canonical variables
for Hamiltonian systems [7, 8]. But the underlying fluid equations are non-
canonical Hamiltonian systems, as is made obvious by the fact that the Euler
equations are highly nonlinear yet their energy function is quadratic [9, 10].
This has motivated the present work, in which we pursue a reformulation of
kinetic wave theory for a number of two-dimensional fluid systems with
quadratic energies based on a particular choice of non-canonical variables.
The practical utility of our choice of variables, which was introduced in a
different context by [11], derives from the existence of a second quadratic
invariant in these systems, which, albeit not sign-definite, greatly
simplifies the wave interaction equations. We leverage these simplifications
into a derivation of scaling laws for the isotropic component of wave spectra
and we present evidence for the importance of these second invariants in
shaping the overall wave spectra in certain situations. Mutatis mutandis, much
of our analysis applies to waves in plasmas as well.

The two-dimensional Boussinesq equations restricted to a vertical $xz$-plane
can be written as

$\displaystyle\Delta\psi_{t}+\left\\{\psi,\Delta\psi\right\\}$
$\displaystyle=-N^{2}\eta_{x}$ (1)
$\displaystyle\eta_{t}+\left\\{\psi,\eta\right\\}$ $\displaystyle=\psi_{x}.$

Here $z$ is the vertical and $x$ is the horizontal coordinate with
corresponding velocities $w$ and $u$, $\psi$ is a stream function such that
$(\psi_{x},\psi_{z})=(w,-u)$ and $-\Delta\psi$ is the vorticity, $\eta$ is the
vertical displacement, $N$ the constant buoyancy frequency and
$\left\\{g,f\right\\}=\partial_{x}g\partial_{z}f-\partial_{z}g\partial_{x}f$.
The vertical buoyancy force $b=-N^{2}\eta$ opposes vertical displacements and
derives from a consideration of potential energy in the presence of gravity
and non-uniform density. It is easily checked that this system has two exact
quadratic invariants: the total energy $E=\int
d\mathbf{x}(-\psi\Delta\psi+N^{2}\eta^{2})$ and the pseudomomentum $P=\int
d\mathbf{x}\,\eta\Delta\psi$. The subtleties associated with the Hamiltonian
point of view of these equations can be appreciated by investigating the
origin of these conservation laws by rewriting (1) as

$\partial_{t}D\phi=\mathcal{J}\frac{\delta E}{\delta\left(D\phi\right)}.$ (2)

Here $\phi^{T}\left(\mathbf{x},t\right)=\left(\psi,\eta\right)$,
$D=\text{diag}(-\Delta,N^{2})$ is a Hermitian semi-positive-definite operator,
and

$\mathcal{J(\phi)}=\frac{1}{2}\begin{pmatrix}\left\\{-\Delta\psi,\cdot\right\\}&\left\\{N^{2}\eta,\cdot\right\\}+N^{2}\partial_{x}\\\
\left\\{N^{2}\eta,\cdot\right\\}+N^{2}\partial_{x}&0\end{pmatrix}$ (3)

is a skew-symmetric operator representing the Poisson structure. This is a
non-canonical Hamiltonian system for the variables $D\phi$ based on the inner
product Hamiltonian function

$E=\left\langle\phi|D\phi\right\rangle=\int\\!d\mathbf{x}\,\phi^{T}D\phi.$ (4)

Energy conservation is then transparently linked to the time translation
symmetry of $\mathcal{J}(\phi)$. The pseudomomentum can be written as

$P=\left\langle D\phi|CD\phi\right\rangle\quad\mbox{with}\quad
C=-\frac{1}{2N^{2}}\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}.$ (5)

Therefore ${\delta P}/{\delta(D\phi)}=2CD\phi$ and

$\mathcal{J}(\phi)\frac{\delta P}{\delta(D\phi)}=-\frac{(D\phi)_{x}}{2},$ (6)

which ensures the invariance of $P$ based on the $x$-translation symmetry of
(4). This suggests interpreting $P$ as a canonical horizontal momentum, even
though it does not agree with the horizontal momentum of the fluid. Actually,
$CD\phi$ is in the kernel of the nonlinear part of $\mathcal{J}$, which
suggests interpreting $P$ as a Casimir of that Poisson degenerate structure
[12, 9]. This degeneracy translates to gauge invariance in terms of Lagrangian
coordinates [13]. Thus, the conservation of $P$ appears Casimir-like based on
the nonlinear dynamics, but momentum-like based on the linear dynamics.
Calling $P$ the ‘pseudomomentum’ is in accordance with established usage in
geophysical fluid dynamics [e.g. §4.3 in 9] and wave–mean interaction theory
[2]. So, whilst the conceptual origins of the two conservation laws for $E$
and $P$ are subtle and subject to interpretation, their actual functional
expressions as quadratic forms $E=\left\langle\phi|D\phi\right\rangle$ and
$P=\left\langle D\phi|CD\phi\right\rangle$ are completely straightforward, and
generalize easily. Following [11], we exploit this by expanding the flow in
variables that diagonalize both $E$ and $P$.

Wave mode expansion. We consider a periodic domain
$\mathbf{x}\in\left[0,L\right]^{2}$ and expand
$\phi(\mathbf{x},t)=\sum_{\alpha}Z_{\alpha}(t)e_{\alpha}(\mathbf{x})$ in terms
of linear wave modes, where $Z_{\alpha}(t)$ are complex scalar wave amplitudes
and the $e_{\alpha}(\mathbf{x})$ are eigenvector functions for the linear part
of (1), i.e.,

$\displaystyle-i\omega_{\alpha}De_{\alpha}$
$\displaystyle={N^{2}}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\partial_{x}e_{\alpha}.$ (7)

The $e_{\alpha}(\mathbf{x})$ are proportional to Fourier modes
$\exp(i\mathbf{k}\cdot\mathbf{x})$ with
$\mathbf{k}=(k_{x},k_{z})\in(2\pi\mathbb{Z}/L)^{2}$. If
$\mathbf{k}=K(\cos\theta,\sin\theta)$ then the dispersion relation is
$\omega=\pm N\cos\theta$. The multi-index $\alpha=(\sigma,\mathbf{k})$
combines branch choice $\sigma=\pm 1$ and wave number $\mathbf{k}$ such that

$\omega_{\alpha=(\sigma,\mathbf{k})}=\sigma N\cos\theta_{k}.$ (8)

With this convention the choice $\sigma=\pm 1$ corresponds to right-going or
left-going waves, respectively, which are vital physical characteristics. The
reality of $\phi$ implies that
$Z_{(\sigma,\mathbf{k})}(t)=Z^{*}_{(\sigma,-\mathbf{k})}(t)$, where the star
denotes complex conjugation. This holds separately within each branch. The
expansion diagonalizes the energy

$\displaystyle
E=\sum_{\alpha}E_{\alpha}=\sum_{\alpha}Z_{\alpha}Z_{\alpha}^{*}$ (9)

and yields the exact equations

$\displaystyle Z_{\alpha,t}+i\omega_{\alpha}Z_{\alpha}$
$\displaystyle=\sum_{\beta,\gamma}\frac{1}{2}V_{\alpha}^{\beta\gamma}Z_{\beta}^{*}Z_{\gamma}^{*}.$
(10)

The interaction coefficients $V_{\alpha}^{\beta\gamma}=\left\langle
e_{\alpha}^{*}|\mathcal{J}(e_{\gamma})e_{\beta}+\mathcal{J}(e_{\beta})e_{\gamma}\right\rangle$
are real, symmetric in upper indices and zero unless
$\mathbf{k}_{\alpha}+\mathbf{k}_{\beta}+\mathbf{k}_{\gamma}=0$. The expansion
also diagonalizes the pseudomomentum

$\displaystyle
P=\sum_{\alpha}P_{\alpha}=\sum_{\alpha}s_{\alpha}Z_{\alpha}Z_{\alpha}^{*},$
(11)

where the horizontal slowness $s_{\alpha}=k_{x}/\omega_{\alpha}$. Hence
$P_{\alpha}=k_{x}/\omega_{\alpha}\,E_{\alpha}$ has the sign of the horizontal
phase (or group) velocity.

The kinetic equation. From (9) and (10) the road to a kinetic equation is
short. The modal wave energy evolves according to

$\dot{E_{\alpha}}=\sum_{\beta,\gamma}V_{\alpha}^{\beta\gamma}\operatorname{Re}\left(Z^{*}_{\alpha}Z^{*}_{\beta}Z^{*}_{\gamma}\right).$
(12)

The kinetic equation describes the evolution of
$n_{\alpha}=\overline{E_{\alpha}}$, where the overbar denotes averaging over a
suitable statistical ensemble. In particular, we average over random Gaussian
initial conditions such that

$\overline{Z_{\beta}^{*}(0)Z_{\alpha}(0)}=\delta_{\alpha\beta}n_{\alpha}(0)$
(13)

is the only nonzero correlation. The standard assumptions and procedural steps
of weak wave turbulence [5, 6] then result in

$\displaystyle\dot{n_{\alpha}}=\pi\\!\\!\int\limits_{\omega_{\alpha\beta\gamma}}\\!\\!V_{\alpha}^{\beta\gamma}\left(V_{\beta}^{\alpha\gamma}n_{\alpha}n_{\gamma}+V_{\gamma}^{\alpha\beta}n_{\beta}n_{\alpha}+V_{\alpha}^{\beta\gamma}n_{\beta}n_{\gamma}\right).$
(14)

Here the joint kinetic limits of big box and long nonlinear times,
$L\rightarrow\infty$ and $t\omega\rightarrow\infty$, were taken. So the
discrete sums in (10) and (12) were replaced by an integral over the resonant
manifold

$\int\limits_{\omega_{\alpha\beta\gamma}}\\!\\!\\!\\!\\!=\\!\int d\beta
d\gamma\,\delta(\omega_{\alpha}+\omega_{\beta}+\omega_{\gamma})\delta\left(\mathbf{k}_{\alpha}+\mathbf{k}_{\beta}+\mathbf{k}_{\gamma}\right),$
(15)

where $\int\\!\\!d\alpha=\sum_{\sigma=\pm 1}\int\\!\\!d\mathbf{k}$. The
kinetic equation (14) is generic for three-wave interactions, but it can be
greatly simplified because of the non-generic additional conservation law for
$P$. The conservation of $E$ and $P$ for every triad (even non-resonant
triads) implies

$\displaystyle
V_{\alpha}^{\beta\gamma}+V_{\beta}^{\alpha\gamma}+V_{\gamma}^{\beta\alpha}$
$\displaystyle=0$ (16) $\displaystyle
s_{\alpha}V_{\alpha}^{\beta\gamma}+s_{\beta}V_{\beta}^{\alpha\gamma}+s_{\gamma}V_{\gamma}^{\beta\alpha}$
$\displaystyle=0,$

respectively [14, 15]. Notably, (16) ensures conservation of $E$ and $P$ for
any projection of (10) onto a truncated set of modes. Viewing (16) as dot
products in $\mathbb{R}^{3}$ means that
$\vec{V}=(V_{\alpha}^{\beta\gamma},V_{\beta}^{\alpha\gamma},V_{\gamma}^{\beta\alpha})$
is orthogonal to both $(1,1,1)$ and $(s_{\alpha},s_{\beta},s_{\gamma})$, which
already determines the direction of $\vec{V}$ uniquely. This makes clear that
any other exact quadratic invariant that is diagonalized by the linear
eigenbasis $\left\\{e_{\alpha}\right\\}$ must be a linear combination of $E$
and $P$. In particular, this implies that the vertical pseudomomentum based on
the other slowness component cannot be exactly conserved by the full dynamics.
However, the kinetic equation is restricted to the resonant manifold
$\omega_{\alpha}+\omega_{\beta}+\omega_{\gamma}=0$, and therefore
$\vec{\omega}=(\omega_{\alpha},\omega_{\beta},\omega_{\gamma})$ is also
orthogonal to both $(1,1,1)$ and $(s_{\alpha},s_{\beta},s_{\gamma})$, the
latter because
$s_{\alpha}\omega_{\alpha}+s_{\beta}\omega_{\beta}+s_{\gamma}\omega_{\gamma}=\hat{\mathbf{x}}\cdot(\mathbf{k}_{\alpha}+\mathbf{k}_{\beta}+\mathbf{k}_{\gamma})=0.$
This means that $\vec{V}$ and $\vec{\omega}$ are parallel to each other, i.e.,

$(V_{\alpha}^{\beta\gamma},V_{\beta}^{\alpha\gamma},V_{\gamma}^{\beta\alpha})=\Gamma_{\alpha\beta\gamma}(\omega_{\alpha},\omega_{\beta},\omega_{\gamma})$
(17)

for some real $\Gamma_{\alpha\beta\gamma}$ totally symmetric in its indices.
For our system

$\displaystyle\Gamma_{\alpha\beta\gamma}=\frac{\left(\sin\theta_{\alpha}+\sin\theta_{\beta}+\sin\theta_{\gamma}\right)}{\sqrt{8}}\left(\sigma_{\alpha}K_{\alpha}+\sigma_{\beta}K_{\beta}+\sigma_{\gamma}K_{\gamma}\right).$
(18)

This changes (14) to

$\dot{n}_{\alpha}={\pi}\\!\\!\int\limits_{\omega_{\alpha\beta\gamma}}\\!\\!\omega_{\alpha}\,\Gamma_{\alpha\beta\gamma}^{2}(\omega_{\alpha}n_{\beta}n_{\gamma}+\omega_{\beta}n_{\alpha}n_{\gamma}+\omega_{\gamma}n_{\alpha}n_{\beta}).$
(19)

Compared to (14), this is a huge simplification and resembles the structure of
kinetic equations derived for canonical Hamiltonian systems (including a
similar entropy function $H(t)=\int d\alpha\log n_{\alpha}$). It is now
apparent that on the resonant manifold additional conservation laws hold
compared to the full system: any component of pseudomomentum is now conserved,
so in $\mathbb{R}^{d}$ there are $d-1$ new conservation laws that are valid
for the kinetic equation. In particular, for the two-dimensional Boussinesq
system the vertical pseudomomentum $P^{z}=\int d\alpha
s^{z}_{\alpha}n_{\alpha}$ is conserved by the kinetic equation but not by the
full flow; here $s^{z}_{\alpha}=k_{z}/\omega_{\alpha}$.

Let us comment on the validity of the kinetic equation with respect to our
initial assumption (13). The off-diagonal correlator
$\overline{Z_{(+,\mathbf{k})}Z_{(-,\mathbf{k})}^{*}}$ has $O(\epsilon^{2})$
fluctuations with frequency $2\omega_{+}(k)$. As long as this beating
frequency does not vanish, the anomalous correlator [16] averages over time
weakly to zero, so that the kinetic equation remains valid. The same holds for
the correlator associated with space homogeneity
$\overline{Z_{(\sigma,\mathbf{k})}^{2}}$. This is discussed in detail in the
supplement together with a detailed derivation of the kinetic equation.

Steady solutions. The frequency (8) and the coefficients (17) are homogeneous
functions of the wavenumbers of degree $w_{\omega}=0$ and $w_{V}=1$,
respectively. That is, for $\lambda>0$,

$\displaystyle\omega_{\left(\sigma,\lambda\mathbf{k}\right)}$
$\displaystyle=\omega_{\left(\sigma,\mathbf{k}\right)}$ (20) $\displaystyle
V_{\lambda\mathbf{k}_{\alpha}}^{\lambda\mathbf{k}_{\beta}\lambda\mathbf{k}_{\gamma}}$
$\displaystyle=\lambda
V_{\mathbf{k}_{\alpha}}^{\mathbf{k}_{\beta}\mathbf{k}_{\gamma}}.$ (21)

Formally at least, we can look for homogeneous solutions to the steady kinetic
equation (19), for which the steady spectrum is then of the separable form

$\displaystyle
n_{\alpha}=n_{\alpha}^{r}(K)n_{\alpha}^{\Omega}\left(\theta_{k}\right).$ (22)

Here $n^{r}_{\alpha}=K^{-w}$ for some suitable $w$. A trivial such formal
steady solution is equipartition of energy such that $w=0$ and
$n^{\Omega}_{\alpha}=$const. This well-known solution has zero spectral flux
of wave energy. Non-zero flux steady solutions are the turbulent solutions,
which are in general not isotropic. However, regardless of their angular
nature, as long as the forcing and dissipation are located at very distinct
scales, an unambiguous spectral energy flux can be defined by integrating (19)
over spheres of radius $K$.

Our result here is that

$n_{\alpha}^{r}=K^{-3}$ (23)

is a steady solution that gives a constant nonzero flux in two distinct cases.
First, $n_{(-,\mathbf{k})}=n_{(+,\mathbf{k})}$, so $P=0$ on average and there
is symmetry between left-going and right-going waves. Second, in cases where
the pseudomomentum is concentrated on just one of the branches, e.g
$n_{(-,\mathbf{k})}=0$, we actually find the two possible power laws

$n^{r}_{E}(K)=K^{-3},\,\,\,n^{r}_{PM}(K)=K^{-3.5},$ (24)

which correspond to energy and pseudomomentum cascades, respectively. In both
cases the isotropic spectrum of two-dimensional internal gravity waves is of
the finite capacity type [17, 18], which shapes the time-dependent self-
similar formation of the steady spectrum and might be of relevance to finite
time singularity formation. See the supplement for a detailed derivation of
these results.

Higher dimensions. We can generalize our results for systems of the form (10)
in $d$-dimensions, with two quadratic invariants: the energy, (9), and the
pseudomomentum, (11). For $P=0$ symmetrically distributed between the
frequency branches, $n_{(\sigma_{i},\mathbf{k})}=n_{(\sigma_{j},\mathbf{k})}$,
we obtain

$n_{E}^{r}(K)=K^{-w^{F}_{E}},\,\,\,\,w^{F}_{E}=w_{V}+d-w_{\omega}/2.$ (25)

Here $w_{V}$ and $w_{\omega}$ are the homogeneity degrees of the interaction
coefficients $V_{\alpha}^{\beta\gamma}$ and of the frequency
$\omega_{\alpha}$. Notably, the isotropic part of the energy spectrum of an
anisotropic system gets an additional contribution of $w_{\omega}/2$ with
respect to the Kolmogorov-Zakharov power law, $w_{KZ}=d+w_{V}-w_{\omega}$, of
Hamiltonian isotropic three-wave interaction systems [5, 19]. Conversely, if
the pseudomomentum is concentrated on one of the branches we find the
additional scaling

$n_{PM}^{r}(K)=K^{-w^{F}_{PM}},\,\,\,\,w^{F}_{PM}=w_{V}+d+(1-2w_{\omega})/2,$
(26)

which corresponds to a cascade of pseudomomentum.

Relevance to other systems. Our results apply to systems described by an
equation of the general form (2) or equivalently (10). Rossby waves in the mid
latitude beta plane, on scales smaller compared to the deformation Radius, are
governed by

$\Delta\psi_{t}+\left\\{\psi,\Delta\psi\right\\}=\beta\partial_{x}\psi$ (27)

where $\psi$ is the stream function on the plane, so $-\Delta\psi$ is the
vorticity. $x$ and $z$ are the zonal and meridional position coordinates and
$f=f_{0}+\beta z$ is the Coriolis parameter. (27) can be written in the form
(2) with
$\phi=\psi,D=-\Delta,\mathcal{J}=\left\\{-\Delta\psi,\cdot\right\\},\mathcal{L}=\beta\partial_{x}$.
It conserves the energy (4) and the pseudomomentum, (5) with $C=1$, which is
the sign-definite enstrophy. Rossby drift waves in plasma are described by a
similar equation. The wave expansion (7) contains only one branch,
$\alpha=\mathbf{k}$, so it is simply the Fourier transform. While the
homogeneity degree of the interaction is the same as for internal gravity
waves, the dispersion relation $\omega_{k}=-\beta\cos\theta/K$ is homogeneous
of degree $w_{\omega}=-1$. Then (25) and (26) give the spectrum

$n_{k}^{r}\propto K^{-3.5},\,\,\,n_{k}^{r}\propto K^{-4.5}$ (28)

for energy and enstrophy cascades, respectively. These scalings agree with the
isotropic part of the unsteady spectra obtained by previous works [20]. Our
work suggests to consider only the isotropic part of these solutions, while
the angular part $n^{\Omega}_{k}$ can be obtained from the kinetic equation
(19) after substituting the radial part $n^{r}_{k}$ from (28).

The dynamics of two-dimensional inertial waves in a vertical plane with
constant Coriolis parameter $f$ is in fact governed precisely by the system
we’ve studied in (2) after the replacement $x\longleftrightarrow
z,f\longleftrightarrow N,\eta\longleftrightarrow v_{2}$. $\psi$ is the stream
function on the vertical plane, $v_{2}$ is the velocity component
perpendicular to the plane. Interestingly, the vertical component of
pseudomomentum is then exactly conserved and is equal to the helicity of the
flow. The isotropic components of energy spectra are given by (23) and (24)
for the limiting cases of zero and sign-definite helicity, respectively.

We are certainly not the first to apply weak wave turbulence approach to study
internal gravity waves [14, 21, 22, 23, 24] and Rossby waves [25, 26, 20, 25].
We believe that our work is the first example of theoretical prediction for
internal gravity waves in $2D$ made by wave turbulence theory and has
experimental practical relevance. Previous studies considered a narrow
spectral range where the homogeneous wave number component is small compared
to the non-homogeneous component. These yield spectral laws with diverging
collision integrals, the divergence of the flux (19), in the case of 3D
internal gravity waves and are not stable in the case of Rossby waves [27] and
hence cannot be physically realized. Our solutions may represent the isotropic
part of a physically relevant solution, while locality (convergence of the
collision integrals) is ensured by obstructions on the complementary angular
part of the energy spectrum $n_{k}^{\Omega}$ (22).

Speculations on inverse energy cascade of unidirectional internal waves. The
decomposition (7) of the field $\phi=\sum_{\alpha}Z_{\alpha}e_{\alpha}$ into
right-going and left-going waves splits the pseudomomentum (11) into positive
and negative components $P=P_{+}+P_{-}$, recall that the positive horizontal
slowness is $s^{x}_{+}=k_{x}/\omega_{+}=K/N$. In the case of waves solely
propagating to the right initially, $Z_{(-,\mathbf{k})}\left(t=0\right)=0$ for
all $\mathbf{k}$, the pseudomomentum is positive at $t=0$, and its time
derivative $\dot{P}_{+}\left(t=0\right)=0$ vanishes as well, as evident from
(12). The persistence of predominantly right-going waves is fixed into the
memory of the system and cannot be forgotten since the pseudomomentum is an
exact invariant; so any generation of leftward propagating waves through
nonlinear interactions must be accompanied by an equal creation of rightward
propagating waves. This has practical relevance for ocean dynamics. For
example, strongly directional internal wave fields arise naturally in the case
of internal tides radiated away from isolated topography structures such as
the Hawaiian ridge [28]. During periods when unidirectional waves dominate,
the system can be effectively described by a kinetic equation truncated to the
positive branch, which shares a similar structure with the kinetic equation of
Rossby waves. While it is well known that the energy of Rossby waves is
transferred, albeit not isotropically, from small to large scales; an inverse
energy cascade of internal gravity waves has never been observed. Ripa [11]
briefly mentions the idea of an inverse energy cascade of internal gravity
waves, adapting the classic dual cascade argument [29]. Similarly to two-
dimensional hydrodynamics [30], the mechanism that drives the inverse energy
cascade of Rossby waves is the existence of a second sign definite quadratic
invariant, the enstrophy, with a density proportional by an isotropic
monotonic function. When unidirectional internal waves dominate, the
pseudomomentum can play a role akin to enstrophy by driving the energy up
scale in the presence of a small or intermediate scale forcing. In this limit
as $P_{+}/P_{-}\rightarrow\infty$, we obtain the scalings for the radial
energy spectra given by (23) and (24) for energy and pseudomomentum cascades,
respectively. For isotropic systems it is widely accepted that the ordering of
the scalings, $3=W^{F}_{E}<W^{F}_{PM}=3.5$, implies an inverse cascade of
energy [31]. While energy transfer in non-isotropic systems can be intricate,
this strongly suggest the existence of inverse energy cascade of internal
waves in the effective description. The extent and time scales of this
approximation for the 2D Boussinesq equation remain a topic for future study.

Adding rotation. The first step towards including both rotation and
stratification in the kinetic equation for internal gravity waves (19) can be
taken whilst retaining the two-dimensional nature of the flow, i.e.,
$\partial_{y}=0$. This involves adding horizontal Coriolis forces to the
momentum equations, which necessitates allowing for a third velocity component
$v_{2}$ in the $y-$direction. The Coriolis forces based on a constant Coriolis
parameter $f$ add only linear terms to the governing equations, so the linear
part of the dynamics is described by the 3x3 operator

$N^{2}\begin{pmatrix}0&0&1\\\ 0&0&0\\\
1&0&0\end{pmatrix}\partial_{x}+f\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}\partial_{z}.$ (29)

The state vector is now $\phi^{T}=(\psi,v_{2},\eta)$ and the Hermitian
diagonal operator is $D=\text{diag}(-\Delta,1,N^{2})$. This generalizes (1) to
include the Coriolis force within the $f$-plane approximation,
$f\hat{z}\times\mathbf{v}$. The energy is
$E=\left\langle\phi|D\phi\right\rangle=\int
d\mathbf{x}\left(-\psi\Delta\psi+v_{2}^{2}+N^{2}\zeta^{2}\right)$. The
dispersion relation is

$\omega_{(\sigma,\mathbf{k})}=\sigma
N\cos\theta_{k}\sqrt{1+f^{2}N^{-2}\tan^{2}\theta_{k}}$ (30)

with $\sigma=0,\pm$. The expansion of $\phi$ in terms of eastward and westward
propagating waves, (7), leads to the kinetic equation (14) with the
interaction coefficients

$\displaystyle V_{\alpha}^{\beta\gamma}$
$\displaystyle=-\frac{\mathbf{k}_{\beta}\times\mathbf{k}_{\gamma}}{2\sqrt{8}K_{\alpha}K_{\beta}K_{\gamma}}\cdot$
(31)
$\displaystyle\left(K_{\beta}^{2}-K_{\gamma}^{2}+f^{2}s_{\alpha}^{z}\left(s_{\gamma}^{z}-s_{\beta}^{z}\right)+N^{2}s_{\alpha}^{x}\left(s_{\gamma}^{x}-s_{\beta}^{x}\right)\right).$

This does not include interactions among and with the zero frequency branch,
also known as the balanced or shear modes. Formally, at the limit of vanishing
frequency, off-diagonal correlators should be added as well to the kinetic
equation. However, weak wave turbulence closure is not expected to remain
valid when shear vortical modes carry the dominant part of the energy [32].
Energy is an exact invariant, so the first constraint in (16) holds. As the
homogeneity degree of the dispersion relation remains $w_{\omega}=0$, our
finding (25), suggests that rotation only modifies the angular component of
the turbulent energy spectrum (22), but leaves the radial component, (23),
unchanged. We note that the conservation of potential vorticity might be used
to simplify the kinetic equation in the presence of both rotation and
stratification. This is studied in a future work.

Conclusion. Our work emphasizes the elegant ramifications in the theoretical
description and in practice of sign-indefinite invariants, which usually do
not get much attention in wave turbulence study. We show that the existence of
a second quadratic invariant, simplifies the kinetic equation and leads to
additional conservation laws on the resonant manifold, which to our knowledge,
were previously unknown in the geophysical community. This simplification
facilitates the derivation of scaling laws for the isotropic component of the
turbulent wave spectra of 2D internal gravity waves and Rossby waves. We show
that there are practical scenarios in which pseudomomentum conservation can
drive an inverse energy cascade of internal gravity waves. On the theoretical
front, our work contributes a different approach to the study of wave
turbulence in non-isotropic systems dominated by three-wave interactions. This
encompasses the application of non-canonical variables for deriving the
kinetic equation, and using variable separation in order to find turbulent
solutions of the kinetic equation.

Acknowledgments. We thank Gregory Falkovich and Vincent Labarre for useful
discussions. This work was supported by the Simons Foundation and the Simons
Collaboration on Wave Turbulence. OB acknowledges additional financial support
under ONR grant N00014-19-1-2407 and NSF grant DMS-2108225. MS acknowledges
additional financial support from the Schmidt Futures Foundation.

## References

  * Andrews  _et al._ [1987] D. G. Andrews, J. R. Holton, and C. B. Leovy,  _Middle atmosphere dynamics_ , 40 (Academic press, 1987).
  * Bühler [2014] O. Bühler,  _Waves and mean flows_ (Cambridge University Press, 2014).
  * Vallis [2017] G. K. Vallis,  _Atmospheric and oceanic fluid dynamics_ (Cambridge University Press, 2017).
  * Whalen  _et al._ [2020] C. B. Whalen, C. De Lavergne, A. C. Naveira Garabato, J. M. Klymak, J. A. MacKinnon, and K. L. Sheen, Internal wave-driven mixing: Governing processes and consequences for climate, Nature Reviews Earth & Environment 1, 606 (2020).
  * Zakharov  _et al._ [2012] V. E. Zakharov, V. S. L’vov, and G. Falkovich,  _Kolmogorov spectra of turbulence I: Wave turbulence_ (Springer Science & Business Media, 2012).
  * Buckmaster  _et al._ [2021] T. Buckmaster, P. Germain, Z. Hani, and J. Shatah, Onset of the wave turbulence description of the longtime behavior of the nonlinear schrödinger equation, Inventiones mathematicae 225, 787 (2021).
  * Zakharov  _et al._ [1987] V. E. Zakharov, A. S. Monin, and L. I. Piterbarg, Hamiltonian description of rossby–blinova baroclinic waves, in  _Doklady Akademii Nauk_ , Vol. 295 (Russian Academy of Sciences, 1987) pp. 1061–1064.
  * Lvov and Tabak [2001a] Y. V. Lvov and E. G. Tabak, Hamiltonian formalism and the garrett-munk spectrum of internal waves in the ocean, Physical review letters 87, 168501 (2001a).
  * Shepherd [1990] T. G. Shepherd, Symmetries, conservation laws, and hamiltonian structure in geophysical fluid dynamics, in  _Advances in Geophysics_ , Vol. 32 (Elsevier, 1990) pp. 287–338.
  * Salmon [1998] R. Salmon,  _Lectures on geophysical fluid dynamics_ (Oxford University Press, 1998).
  * Ripa [1981] P. Ripa, On the theory of nonlinear wave-wave interactions among geophysical waves, Journal of Fluid Mechanics 103, 87 (1981).
  * Olver [1982] P. J. Olver, A nonlinear hamiltonian structure for the euler equations, Journal of Mathematical Analysis and Applications 89, 233 (1982).
  * Olver [1993] P. J. Olver,  _Applications of Lie groups to differential equations_ , Vol. 107 (Springer Science & Business Media, 1993).
  * Hasselmann [1966] K. Hasselmann, Feynman diagrams and interaction rules of wave-wave scattering processes, Reviews of Geophysics 4, 1 (1966).
  * Vanneste [2005] J. Vanneste, Wave interactions (Springer Vienna, Vienna, 2005) pp. 69–94.
  * Zaleski  _et al._ [2020] J. Zaleski, M. Onorato, and Y. V. Lvov, Anomalous correlators in nonlinear dispersive wave systems, Physical Review X 10, 021043 (2020).
  * Connaughton and Newell [2010] C. Connaughton and A. C. Newell, Dynamical scaling and the finite-capacity anomaly in three-wave turbulence, Physical Review E 81, 036303 (2010).
  * Krstulovic and Nazarenko [2023] G. Krstulovic and S. Nazarenko, Initial evolution of 3d turbulence en route to the kolmogorov state: emergence and transformations of coherent structures, self-similarity and instabilities, arXiv preprint arXiv:2310.01121 (2023).
  * Zakharov [1965] V. E. Zakharov, Weak turbulence in media with a decay spectrum, Journal of Applied Mechanics and Technical Physics 6, 22 (1965).
  * Pelinovsky [1978] E. Pelinovsky, Wave turbulence on beta-plane, Okeanologiya 18, 192 (1978).
  * McComas and Bretherton [1977] C. H. McComas and F. P. Bretherton, Resonant interaction of oceanic internal waves, Journal of Geophysical Research 82, 1397 (1977).
  * Müller  _et al._ [1986] P. Müller, G. Holloway, F. Henyey, and N. Pomphrey, Nonlinear interactions among internal gravity waves, Reviews of Geophysics 24, 493 (1986).
  * Caillol and Zeitlin [2000] P. Caillol and V. Zeitlin, Kinetic equations and stationary energy spectra of weakly nonlinear internal gravity waves, Dynamics of atmospheres and oceans 32, 81 (2000).
  * Lvov and Tabak [2001b] Y. V. Lvov and E. G. Tabak, Hamiltonian formalism and the garrett-munk spectrum of internal waves in the ocean, Physical review letters 87, 168501 (2001b).
  * Monin and Piterbarg [1987] A. S. Monin and L. I. Piterbarg, The kinetic equation for rossby–blinova waves, in  _Doklady Akademii Nauk_ , Vol. 295 (Russian Academy of Sciences, 1987) pp. 816–820.
  * Connaughton  _et al._ [2015] C. Connaughton, S. Nazarenko, and B. Quinn, Rossby and drift wave turbulence and zonal flows: the charney–hasegawa–mima model and its extensions, Physics Reports 604, 1 (2015).
  * Balk and Nazarenko [1990] A. Balk and S. Nazarenko, On the physical realizability of anisotropic kolmogorov spectra of weak turbulence, Sov. Phys. JETP 70 (1990).
  * Smith and Young [2003] S. G. L. Smith and W. Young, Tidal conversion at a very steep ridge, Journal of Fluid Mechanics 495, 175 (2003).
  * Fjørtoft [1953] R. Fjørtoft, On the changes in the spectral distribution of kinetic energy for twodimensional, nondivergent flow, Tellus 5, 225 (1953).
  * Kraichnan [1967] R. H. Kraichnan, Inertial ranges in two-dimensional turbulence, The Physics of Fluids 10, 1417 (1967).
  * Nazarenko [2011] S. Nazarenko,  _Wave turbulence_ , Vol. 825 (Springer Science & Business Media, 2011).
  * Labarre  _et al._ [2023] V. Labarre, P. Augier, G. Krstulovic, and S. Nazarenko, Internal gravity waves in stratified flows with and without vortical modes, arXiv preprint arXiv:2303.01570 (2023).