#  Non - relativistic string in Newton - Cartan background

Sk. Moinuddin $\,{}^{\rm a,b}$, Pradip Mukherjee $\,{}^{\rm a,c}$

aDepartment of Physics, Barasat Government College,Barasat, India

<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
##  1 Introduction

Gravitation is the weakest of the fundamental forces of nature.
Notwithstanding this, the physics in the large scale is overwhelmingly
influenced by the gravitational interaction. The reason is, gravitation is
always attractive . Einstein in his celebrated General (Theory of ) relativity
(GR) formulated the gravitational interaction as the change of curvature of
spacetime . As a result free particles follow the shortest path (geodesic) in
this curved manifold. So Einstein’s equations that determine the curvature
produced by a certain mass energy distribution. The particle (or string) under
study are the test particles. The test particles are so chosen that they do
not disturb the geometry in any way [1]. This framework should not be violated
in any limiting procedure. The non relativistic theory which we are
considering should not be an exception.

Physical theories are based on space and time. In Newtonian concept time flows
from past to future in an absolute sense and the physical space is relative.
The geometry of space is euclidean and one can set up a Cartesian space system
in which the space metric is $\delta_{(}ij)$. This is the non relativistic
flat space. Gravitation of massive bodies produce curvature [1] in the free
space. The similarity and differences with relativity theory are obvious.
Also, the more accurate theory of special relativity in the non relativistic
limit $(v<<c)$ should agree with non relativistic results. The boost invariant
non relativistic spacetime should also follow as a limit of the relativistic
Minkowski spacetime. Indeed, one can demonstrate the transition by examining
the structure of the metric in the $c\to\infty$ limit to pass into a
collection of singular metrics of rank 1 and 3 respectively [26]. Different
non relativistic geometries may follow from the limiting procedure.

Both classical theory of gravitation and Einstein’s general relativity (GR)
obey the principle of equivalence by dint of which gravitational effects may
be annulled locally in the $\epsilon$ \- neighborhood of a point. This was
exploited by Cartan to discover the geometric formulation of Newtonian
gravity. The corresponding space time manifold is known as Newton - Cartan
spacetime [1]. Research on the Newton Cartan spacetime is continued over a
long time [9]. However, these theories are mostly formulated as second order
metric theory, following GR. The vielbein (tetrad ) formulation of the NC
gravity was less emphasized. Some authors did use the vielbein structures but
an out and out discussion in terms of the vielbeins and spin connections is
scarcely found in the literature.

The scenario has changed in last couple of decades when non relativistic field
theories were found to be important in the context of strongly correlated
electrons in mesoscopic physics and soft condensed matter systems [24]. The
requirements are now coupling different particle and field theories to non
relativistic curved space time , This coupling is very simple for the
relativistic theories. One has to replace the ordinary derivatives by the
corresponding co-variant derivatives and replace the measure of integration by
the invariant measure . Now substituting the Minkowski metric $\eta_{\mu\nu}$
by the metric $g_{\mu\nu}$ of the curved spacetime. In case of the non
relativistic geometry the metric is degenerate. There is a spatial metric
$h^{\mu\nu}$ of rank three and a one form $\tau_{\mu}$. Since the metric
structure is singular they can not be inverted. Thus the theory formulated by
vielbeins become important. Now in GR, the metric can be factored into the
vielbeins. This is not available in the non relativistic counterpart. Rather,
the formulation in the non relativistic spacetime is far easier in the first
order approach with the geometry given by the tetrad and the spin connections
, as their introduction do not depend upon the metric. So the question is how
to formulate the first order theory, without taking recourse of the metric?

Consider a theory invariant under the Poincare transformations in the
Minkowski space. An interesting correspondence is unveiled by localizing
symmetry transformation parameters of a theory under the Poincare group of
transformations in the Minkowski space. Fields would transform formally in the
same way but the derivatives transform differently as the parameters now are
no more constants. The symmetry is regained when appropriate gauge fields are
included which transform appropriately. These transformation of the fields
have one to one correspondence with the vielbeins and spin connections of the
Riemann Cartan manifold. The theory thus emerges is the famous Poincare gauge
theory (PGT) 111Of course the target of PGT was to provide a quantum theory of
gravity which is so far unachieved.. Naturally, one would like to apply the
idea of PGT to the non relativistic arena. In the non relativistic problem, as
we have pointed out already, the metric is degenerate. The gauging of symmetry
approach seems to be an elegant and robust method for coupling different
matter theories with the curved back ground produced by gravity. This was
first demonstrated for the nonrelativistic Schroedinger field in [10] which
resolved the paradoxes in the treatment of the same problem elsewhere [13].
Soon the connection of the localisation process with the emergence of the
Newton - Cartan geometry was shown [12] and extension of the algorithm to
include external electromagnetic interaction was presented [11, 10]. Finally
the gauging of the symmetry approach was formalised as Galilean gauge theory
(GGT) [10]. Note that the Galilean symmetry is inbuilt in the algorithm of
GGT. This new theory is suitable for the purpose due to the following
features:

  1. 1.

The algorithm developed in the theory can be useful for any theory which is
symmetric in flat Galilean spacetime. The method of approach automatically
carries the symmetry along with it. Thus failure to reproduce the flat theory
in the appropriate limit , so common in other approaches, has no place in this
approach.

  2. 2.

The whole calculations are done by a set of rules, derived once for all .
There is no fine adjustments during the calculation of a particular problem.

  3. 3.

The spacetime emerging from our analysis is the Newton- Cartan spacetime. So
far we have considered examples from field theory, particle models. In this
paper we will provide our results for string model.

##  2 Action for Non relativistic Particle in Curved Background

We start with a short review of the calculations and results for the
nonrelativistic particle model treated by the Galilean gauge theory (GGT)
[33]. The parametrized action for a non relativistic particle in $3$
dimensional Euclidean space and absolute time is given by,

$S=\int\dfrac{1}{2}m\dfrac{\dfrac{dx^{a}}{d\lambda}\dfrac{dx^{a}}{d\lambda}}{\bigg{(}\dfrac{dx^{0}}{d\lambda}\bigg{)}}~{}d\lambda$
(1)

where the index $a$ denotes a space index. The action is invariant under the
global Galilean transformations,

$x^{\mu}\to
x^{\mu}+\xi^{\mu};\xi^{0}=-\epsilon,\xi^{k}=\eta^{k}-v^{k}t;\eta^{k}=\omega^{k}{}_{l}x^{l}+\epsilon^{k}$
(2)

This invariance is ensured by the transformations

$\delta\dfrac{dx^{0}}{d\lambda}=\dfrac{d}{d\lambda}(\delta
x^{0})=-\dfrac{d\epsilon}{d\lambda}=0$ (3)

as $\epsilon$ is constant and,

$\delta\dfrac{dx^{k}}{d\lambda}=w^{k}_{j}\dfrac{dx^{j}}{d\lambda}-v^{k}\dfrac{dx^{0}}{d\lambda}$
(4)

which can be checked easily. The Lagrangian (1) changes by,

$\delta L=-\dfrac{d}{d\lambda}\bigg{(}mv^{k}\dfrac{dx^{k}}{d\lambda}\bigg{)}$
(5)

due to (3) and (4). The change of the action (1) is then a boundary terms
only, see (5). The same equations of motion follow from both the original and
the transformed action. So the theory is invariant under the global Galilean
transformations.

To localize the symmetry of the action (1) according to GGT
,$\frac{dx^{\alpha}}{d\lambda}$ is now subsituted by the covariant derivatives
$\frac{Dx^{\alpha}}{d\lambda}$, where

$\dfrac{Dx^{\alpha}}{d\lambda}=\dfrac{dx^{\nu}}{d\lambda}\Lambda^{\beta}_{\nu}\partial_{\beta}x^{\alpha}=\dfrac{dx^{\nu}}{d\lambda}\Lambda^{\alpha}_{\nu}$
(6)

Here $\Lambda^{\beta}_{\nu}$ are a set of new compensating (gauge) fields.
This covariant derivative was elaborately derived. We are not repeating it
here. Instead let us stress that this will act as a building block for the
string, as the latter is viewed here as an extension of the covariant particle
model in [33], where it was derived first in this approach. Note that the
localisation procedure can be tuned smoothly to restore global Galilean
symmetry. In this limit the covariant derivatives
$\frac{Dx^{\alpha}}{d\lambda}$ must then go to the ordinary derivatives
$\frac{dx^{\mu}}{d\lambda}$. This sets the following condition

$\Lambda^{\alpha}{}_{\mu}\longrightarrow\delta^{\alpha}_{\mu}$ (7)

This at once shows that $\Lambda^{\alpha}{}_{\mu}$ is non singular i.e. the
corresponding matrix is invertible. We will denote the inverse by
$\Sigma^{\mu}_{\alpha}$. We will further see that this observation is
instrumental in the geometric interpretation of our theory.

Now the transformations of the new gauge fields should ensure that the
’covariant derivatives’ will transform under the local Galilean transfrmations
in the same form as the usual derivatives do under the global Galilean
transformations. Then the new theory obtained by replacing the ordinary
derivatives by the ‘covariant derivatives’ will be invariant under the local
Galilean transformations. This is the essence of GGT [10, 12]. Using (3) and
(4) the transformation of the covariant derivatives follows,

$\delta\dfrac{Dx^{0}}{d\lambda}=0$ (8)

and, likewise for the space part,

$\delta\dfrac{Dx^{k}}{d\lambda}=w^{k}_{j}\dfrac{Dx^{j}}{d\lambda}-v^{k}\dfrac{Dx^{0}}{d\lambda}$
(9)

Exploiting these relations the transformations of the newly introduced fields
are completely specified. They are given by [12],

$\displaystyle\delta\Lambda_{0}^{0}$ $\displaystyle=$
$\displaystyle\dot{\epsilon}\Lambda_{0}^{0}$
$\displaystyle\delta\Lambda_{i}^{a}$ $\displaystyle=$
$\displaystyle\omega_{b}^{a}\Lambda_{i}^{b}-\partial_{i}\xi^{k}\Lambda_{k}^{a}$
$\displaystyle\delta\Lambda_{0}^{a}$ $\displaystyle=$
$\displaystyle\dot{\epsilon}\Lambda_{0}^{a}-v^{a}\Lambda_{0}^{0}-\partial_{0}\xi^{k}\Lambda_{k}^{a}+\omega_{b}^{a}\Lambda_{0}^{b}$
(10)

while the remaining field $\Lambda_{i}^{0}$ simply vanishes.

Hence the cherished action is given by,

$S=\int\dfrac{1}{2}m\dfrac{\dfrac{Dx^{a}}{d\lambda}\dfrac{Dx^{a}}{d\lambda}}{\bigg{(}\dfrac{Dx^{0}}{d\lambda}\bigg{)}}~{}d\lambda$
(11)

obtained from (1), substituting $\frac{dx^{a}}{d\lambda}$ by
$\frac{Dx^{a}}{d\lambda}$. The above derivation is the standard procedure of
GGT which leads to the theory (11) invariant under transformations formally
similar to (2). But the parameters of transformations are not constants but
vary in a special manner with space and time, completely in unison with the
privileged role of time in Newtonian theory.

##  3  Non-relativistic Nambu - Goto action for the bosonic string

The string is an extension of the particle model, Both are relativistic
objects [33] and a particular type of nonrelativistic limit is to be taken to
obtain the coupling of the nonrelativistic theories with gravity. The
algorithm of GGT is remarkable successful in formulating the coupling in
different particle models [10]. So one would expect that the same algorithm
will be useful for coupling the strings though we must be very careful on
account of some very important differences between the two models. Unlike the
particle which is represented by a point, the string is an one dimensional
object which is described by a parameter $\sigma$. So during its evolution it
traces a two dimensional world sheet. This surface defines a two dimensional
foliation of space time. It is mapped by two coordinates, $\tau$ and $\sigma$,
where $\tau$ is time-like and $\sigma$ is space-like. As the string is
essentially relativistic, the low energy excitation do not affect the string
world sheet The relativistic Nembu Goto action of the bosonic string is given
by,

$S_{\rm{NG}}=-N\int{d\sigma}{d\tau}\sqrt{-\det h_{ij}}$ (12)

where $h_{ij}$ is the metric induced by the target space and given by,

$h_{ij}=\eta_{\mu\nu}\partial_{i}X^{\mu}\partial_{j}X^{\nu}$ (13)

where $X^{\mu}=X^{\mu}(\tau,\sigma)$ and
$\eta_{\mu\nu}={\rm{diag}1,\quad-1,\quad-1....}$ is the Lorentzian metric in
the target space. String (12) is a relativistic object but if we consider low
energy phenomenology, then the effects in the target space are non
relativistic. Let the dimension of the embedding space be (D+1). At a given
time the string intersects the embedding space along a line. We take this line
as the $X^{1}$ coordinate line. Then $X^{0}$ and $X^{1}$ are longitudinal to
the string and the rest is transverse. Note that we could take any of
$X^{2}$,……. $X^{D}$ in place of $X^{1}$. So no particular gauge choice is
associated with this. By taking $c\to\infty$ limit, we get the non
relativistic Lagrangian,

$\displaystyle{\mathcal{L}}_{NG}=-N{\left(2\epsilon_{\mu\nu}\dot{X}^{\mu}{X^{\prime}}{}^{\nu}\right)}^{-1}\left(\dot{X}^{\mu}X^{a\prime}-\dot{X}^{a}X^{\mu\prime}\right)^{2}$
(14)

and the action is given by

$S_{NG}=\int d\sigma d\tau{\mathcal{L}}_{NG}$ (15)

with ${\mathcal{L}}_{NG}$ given by (14) is Galilean invariant. The non
relativistic lagrangian (14) can be written as

$\displaystyle{\mathcal{L}}_{NG}=-N{\left(\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\right)}^{-1}\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{a}}{\partial\sigma_{\beta}}\right)^{2}$ (16)

Where $\alpha,\beta=0,1$ and $\sigma_{1}=\tau$ and $\sigma_{2}=\sigma$.

With this we finish the review of the method applied to a generally covariant
particle model and introduce the main results of our construction of the
string (bosonic) action following the particle action, we are in a position to
introduce gravity. The string being an essentially relativistic object, the
low energy excitations of the string are confined in the direction transverse
to the string world sheet. The curvature produced by gravity are only to be
considered. We describe the procedure to non - relativistic.

##  4  Non-relativistic bosonic string in curved background from galilean
gauge theory

Coupling a Non relativistic string theory with gravity has been found to be a
difficult task because of the peculiar geometry of the non- relativistic
string embedded in curved space time. The string is inherently a relativistic
object. So the 2 dimensional slice of the space time produced by the motion of
the sting (the world sheet) has a Minkowski metric structure where as the
transverse bulk has an Euclidean structure. The point is that the non
relativistic description for the string is relevant for low energy
excitations, thus in interaction with gravity the world sheet is not affected.
Remember that in formulating the non relativistic action for flat space time
we have enforced the condition $\omega^{1a}=0$ , where $\omega$ is the spatial
rotation parameter , the coordinate axis $X^{1}$ is longitudinal and $X^{a}$
are transverse to world sheet .

Once the peculiarity of the non relativistic string geometry is understood it
is simple to write down the corresponding action for such string coupled with
curved manifold , thanks to the algorithm of GGT . We will thus replace the
ordinary derivative by covariant derivative. Note that this replacement is
with in the transverse part of the manifold.

Hence the lagrangian of non relativistic bosonic string (16) in curved
background is given by

$\displaystyle{\mathcal{L}}_{NG}=-N{\left(\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\right)}^{-1}\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\right)^{2}$
(17)

Which is a direct generalization of the nonrelativistic particle model [33].
Explicitly

$\displaystyle\dfrac{DX^{a}}{d\sigma_{\beta}}$ $\displaystyle=$
$\displaystyle\dfrac{dX^{\rho}}{d\sigma_{\beta}}\Lambda^{a}_{\rho}$
$\displaystyle\delta\dfrac{DX^{a}}{d\sigma_{\beta}}$ $\displaystyle=$
$\displaystyle\dfrac{DX^{j}}{d\sigma_{\beta}}\omega^{a}_{j}-\dfrac{DX^{0}}{d\sigma_{\beta}}u^{a}$
$\displaystyle\delta\dfrac{\partial X^{\mu}}{\partial\sigma_{\alpha}}$
$\displaystyle=$ $\displaystyle\dfrac{\partial
X^{\nu}}{d\sigma_{\alpha}}\omega^{\mu}_{\nu}$ (18)

Note that those degrees of freedom transverse to the string world sheet is
interacted by the nonrelativistic gravity. We have already discussed this
issue in the above.

##  5 General covariance of the theory

Our next task is to show that the (17) is invariant under local Galilean
transformation. This is equivalent to that the diffeomorphism invariance

$X^{\mu}\to X^{\mu}+\xi^{\mu}$ (19)

where

$\displaystyle\xi^{\lambda}$ $\displaystyle=$
$\displaystyle\epsilon^{\lambda}+w^{\lambda}_{\omega}x^{\omega}~{}~{},~{}~{}\lambda,\omega=0,1$
$\displaystyle\xi^{a}$ $\displaystyle=$
$\displaystyle\epsilon^{a}+w^{a}_{j}x^{j}-u^{a}t~{}~{},~{}~{}a=2........D.$
(20)

holds.

Since the transformation parameters are abitrary functions of coordinates and
time the diffeomorphism (19) is totally arbitrary . The variation of
Lagrangian due to the spacetime diffeomorphism is ,

$\displaystyle\delta{\mathcal{L}}_{NG}$ $\displaystyle=$
$\displaystyle-N{\left(\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\right)}^{-1}\bigg{[}\delta\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\right)^{2}\bigg{]}$
(21)
$\displaystyle\qquad\qquad-N\bigg{[}\delta{\left(\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\right)}^{-1}\bigg{]}\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\right)^{2}$

Now the calculation of (21) is very simple but quite lengthy , so we calculate
it in term by term. The 2nd term of the (21)

$\displaystyle\delta{\left(\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\right)}^{-1}$ $\displaystyle=$
$\displaystyle-\frac{\epsilon_{\mu\nu}\sigma_{\alpha\beta}\bigg{[}\delta\bigg{(}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\bigg{)}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}+\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\delta\bigg{(}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\bigg{)}\bigg{]}}{\bigg{(}\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\bigg{)}^{2}}$
$\displaystyle\qquad\qquad=-\frac{\epsilon_{\mu\nu}\sigma_{\alpha\beta}\bigg{[}\omega^{\mu}_{\nu}\frac{\partial
X^{\nu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}+\omega^{\nu}_{\mu}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\mu}}{\partial\sigma_{\beta}}\bigg{]}}{\bigg{(}\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\bigg{)}^{2}}$

as $\epsilon_{\mu\nu}$ is anti symmetric thus

$\displaystyle\delta{\left(\epsilon_{\mu\nu}\sigma_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma_{\beta}}\right)}^{-1}=0$ (22)

now 1st term of the (21)

$\displaystyle\delta\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\right)^{2}$
$\displaystyle=$ $\displaystyle\bigg{(}2\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\bigg{)}\bigg{[}\delta\bigg{(}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\bigg{)}\frac{DX^{a}}{d\sigma_{\beta}}\epsilon_{\alpha\beta}+\delta\bigg{(}\frac{DX^{a}}{d\sigma_{\beta}}\bigg{)}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\epsilon_{\alpha\beta}\bigg{]}$
$\displaystyle\qquad\qquad=\bigg{(}2\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\bigg{)}\bigg{[}\bigg{(}\frac{\partial
X^{\nu}}{\partial\sigma_{\alpha}}\omega^{\mu}_{\nu}\bigg{)}\frac{DX^{a}}{d\sigma_{\beta}}\epsilon_{\alpha\beta}+$
$\displaystyle\qquad\qquad\bigg{(}\dfrac{DX^{j}}{d\sigma_{\beta}}\omega^{a}_{j}-\dfrac{DX^{0}}{d\sigma_{\beta}}u^{a}\bigg{)}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\epsilon_{\alpha\beta}\bigg{]}$

putting the value of $\dfrac{DX^{j}}{d\sigma_{\beta}}$ and
$\dfrac{DX^{0}}{d\sigma_{\beta}}$ we get

$\displaystyle\delta\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\right)^{2}$
$\displaystyle=$ $\displaystyle\bigg{(}2\epsilon_{\lambda\rho}\frac{\partial
X^{\mu}}{\partial\sigma_{\lambda}}\frac{dX^{m}}{d\sigma_{\rho}}\Lambda^{a}_{m}\bigg{)}\bigg{[}\epsilon_{\alpha\beta}\omega^{\mu}_{\nu}\frac{\partial
X^{\nu}}{\partial\sigma_{\alpha}}\frac{dX^{i}}{d\sigma_{\beta}}\Lambda^{a}_{i}+$
(23)
$\displaystyle\qquad\qquad\bigg{(}\omega^{a}_{j}\dfrac{dX^{l}}{d\sigma_{\beta}}\Lambda^{j}_{l}-u^{a}\dfrac{dX^{n}}{d\sigma_{\beta}}\Lambda^{o}_{n}\bigg{)}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\epsilon_{\alpha\beta}\bigg{]}$

Now $\Lambda^{o}_{n}=0$ and so the equation (23) becomes

$\displaystyle\delta\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\right)^{2}$
$\displaystyle=$ $\displaystyle\bigg{(}2\epsilon_{\lambda\rho}\frac{\partial
X^{\mu}}{\partial\sigma_{\lambda}}\frac{dX^{m}}{d\sigma_{\rho}}\Lambda^{a}_{m}\bigg{)}\bigg{[}\epsilon_{\alpha\beta}\omega^{\mu}_{\nu}\frac{\partial
X^{\nu}}{\partial\sigma_{\alpha}}\frac{dX^{i}}{d\sigma_{\beta}}\Lambda^{a}_{i}+\epsilon_{\alpha\beta}\omega^{a}_{j}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\dfrac{dX^{l}}{d\sigma_{\beta}}\Lambda^{j}_{l}\bigg{]}$
(24)

the 1st term of (24)

$\displaystyle
2\epsilon_{\lambda\rho}\epsilon_{\alpha\beta}\omega^{\mu}_{\nu}\frac{\partial
X^{\mu}}{\partial\sigma_{\lambda}}\frac{\partial
X^{\nu}}{\partial\sigma_{\alpha}}\frac{dX^{m}}{d\sigma_{\rho}}\frac{dX^{i}}{d\sigma_{\beta}}\Lambda^{a}_{m}\Lambda^{a}_{i}$
$\displaystyle=$
$\displaystyle\bigg{(}\epsilon_{\lambda\rho}\epsilon_{\alpha\beta}\omega^{\mu}_{\nu}\frac{\partial
X^{\mu}}{\partial\sigma_{\lambda}}\frac{\partial
X^{\nu}}{\partial\sigma_{\alpha}}\frac{dX^{m}}{d\sigma_{\rho}}\frac{dX^{i}}{d\sigma_{\beta}}\Lambda^{a}_{m}\Lambda^{a}_{i}+$
$\displaystyle\qquad\qquad\epsilon_{\lambda\rho}\epsilon_{\alpha\beta}\omega^{\mu}_{\nu}\frac{\partial
X^{\mu}}{\partial\sigma_{\lambda}}\frac{\partial
X^{\nu}}{\partial\sigma_{\alpha}}\frac{dX^{m}}{d\sigma_{\rho}}\frac{dX^{i}}{d\sigma_{\beta}}\Lambda^{a}_{m}\Lambda^{a}_{i}\bigg{)}$

swapping indices and using the different symmetry property we find that it
cancels the second term of the same. Thus we see that the first term of (24)
vanish. Similarly we can proved that the second terms of (24) also vanish.
Thus

$\displaystyle\delta\left(\epsilon_{\alpha\beta}\frac{\partial
X^{\mu}}{\partial\sigma_{\alpha}}\frac{DX^{a}}{d\sigma_{\beta}}\right)^{2}=0$
(25)

So from (25) and (22) we see that

$\displaystyle\delta{\mathcal{L}}_{NG}=0$ (26)

Thus the curved background Lagrangian (17) is invariant under diffeomorphism
due to local Galilean transformations. The action (17) is quite satisfactory.
The flat limit poses no problems. In this limit the vielbeins reduce to the
Kronecker deltas and using (18) we can easily show that The action (17)
reduces to the NR string action (16) .

##  6 Connection with geometry

In this section we will elaborate how the above action belongs to the Newton
Cartan geometry. This is important because several assumptions were made
regarding the geometry of the string. We will remember that a follow up
assumption was the excitations of the string are entirely confined to the
brane transverse to the string world sheet. Application of the Galilean gauge
theory further complicates the issue. A look at equation (19) is enough to
rouse our scepticism. However, this scepticism is quenched very soon when we
proved that the flat model [36] is invariant under Galilean transformations.
It is intuitively obvious that the curved space with which the string model is
coupled is the Newton Cartan spacetime. Still, to appreciate the result the
following remarks may be useful. The strategy is this. We will give the
principal results for the localization of a generic model endowed with the all
the symmetries of the Galilean group 222but not any central charge. The
resulting curved space time is demonstrated analytically to be the N - C
spacetime [10]. Theories which have Galilean symmetry in the flat limit when
gauged may have fewer transformations but those transformations form an
illustration of the Newton Cartan algebra.

The GGT is primarily a gauge theory in flat space time . Now to account for
the local degrees of freedom a local coordinate basis is set up at every point
where the basis vectors trivially parallel to the respective global basis
vectors. The transformations of the gauge fields introduced during
localization have a suggestive form which are begging for a geometrical
interpretation [28] . The transformations obtained here (10) are the same as
the corresponding one in the general results. These transformations are just
the same as the vielbeins and spin connections of a curved manifold. The new
fields ${\Lambda_{\mu}}^{\alpha}$ may then be reinterpreted as inverse
vielbeins in a general manifold charted by the coordinates $x^{\mu}$
connecting the (global) spacetime coordinates with the local coordinates. A
crucial requirement for this is the existence of an inverse which is met above
by the continuity of the local to global invariance (see the equation (7) and
the discussion around it ). 333A word about our notation: indices from the
beginning of the alphabet ($\alpha,\beta$ etc. or $a,b$ etc) denote the local
basis while those from the middle ($\mu,\nu$ etc. or $i,j$ etc.) indicate the
global basis; Greek indices denote space-time while only space is given by the
Latin ones. Repeated indices denote a summation. The geometric
reinterpretation is then seen from the above transformations see, for example
(10) where the local indices (denoted by $a$) are Lorentz rotated, while the
global indices (denoted by $i$) are coordinate transformed.

The generators of the space time translations are ${\Lambda_{\mu}}^{\alpha}$
and those of rotations are ${B^{\mu}}_{ab}$ and boosts are $B^{\mu}_{a0}$ . We
observe that the rotation generators have the usual anti symmetry.

We introduce ${\Sigma_{\alpha}}^{\mu}$ as the inverse of
$\Lambda_{\mu}{}^{\beta}$ . So

${\Sigma_{\alpha}}^{\mu}\Lambda_{\mu}{}^{\beta}=\delta^{\beta}_{\alpha},~{}~{}~{}{\Sigma_{\alpha}}^{\mu}\Lambda_{\nu}{}^{\alpha}=\delta^{\mu}_{\nu}$
(27)

The spatial part $\Lambda_{k}{}^{a}$ is the inverse of $\Sigma_{a}{}^{k}$ as
may be seen in (27). Note that we are denoting the local coordinates by the
initial Greek letters i.e. $\alpha,\beta$ etc. whereas the global coordinates
are denoted by letters from the middle of the Greek alphabet i.e. $\mu,\nu$
etc. From the definitions of ${\Sigma_{\alpha}}^{\mu}$ and
$\Lambda_{\mu}{}^{\alpha}$ and the transformations of the various fields
involved we can obtain the corresponding transformation laws. Thus,

$\displaystyle\delta_{0}{\Sigma_{0}}^{k}$
$\displaystyle=-\xi^{\nu}{\partial_{\nu}{\Sigma}_{0}}^{k}+{\Sigma_{0}}^{\nu}\partial_{\nu}\xi^{k}-v^{b}{\Sigma_{b}}^{k}$
$\displaystyle\delta_{0}{\Sigma_{a}}^{k}$
$\displaystyle=-\xi^{\nu}{\partial_{\nu}{\Sigma}_{a}}^{k}+{\Sigma_{a}}^{\nu}\partial_{\nu}\xi^{k}-\lambda_{a}{}^{b}{\Sigma_{b}}^{k}$
(28)

It has been proved that the 4-dim space time obtained in this way above is the
Newton-Cartan manifold. This is done by showing that the metric formulation of
our theory contains the same structures and satisfy the same structural
relations as in NC spacetime [12]. the invariance of the theory is retained
after localization if the newly introduced basic fields transform as

$\displaystyle\delta_{0}\Sigma_{0}{}^{0}$ $\displaystyle=$
$\displaystyle-\zeta^{\nu}\partial_{\nu}\Sigma_{0}{}^{0}+\Sigma_{0}{}^{\nu}\partial_{\nu}\zeta^{0}$
$\displaystyle\delta_{0}\Sigma_{0}{}^{k}$ $\displaystyle=$
$\displaystyle-\zeta^{\nu}\partial_{\nu}\Sigma_{0}{}^{k}+\Sigma_{0}{}^{\nu}\partial_{\nu}\zeta^{k}+u^{a}\Sigma_{a}{}^{k}$
$\displaystyle\delta_{0}\Sigma_{a}{}^{k}$ $\displaystyle=$
$\displaystyle-\zeta^{\nu}\partial_{\nu}\Sigma_{a}{}^{k}+\Sigma_{a}{}^{\nu}\partial_{\nu}\zeta^{k}-\omega_{a}{}^{b}\Sigma_{b}{}^{k}$
$\displaystyle\delta_{0}B_{\mu}$ $\displaystyle=$
$\displaystyle-\zeta^{\nu}\partial_{\nu}B_{\mu}-\partial_{\mu}\zeta^{\nu}B_{\nu}-\frac{1}{2}\partial_{\mu}\omega^{ab}\sigma_{ab}+\partial_{\mu}u^{a}{x_{a}}$
(29)

Equation (29) give the complete structure of the Newton - Cartan space time.
What we would like to emphasise is in a particular example not all the fields
will appear. Which fields would be included on the dynamics of the particular
model in question. The examples of the non - relativistic particle and bosonic
string are quiet instructive. Both carry the sense of time along the path. The
basic dynamical variables do not carry any intrinsic direction attached with
it. The required transformations are now [33]

$\epsilon\to\epsilon(t)\,\,\,;\,\,\,\epsilon^{k},\,\,\,w^{k}_{j},\,\,\,u^{k}\to
f(x,t)$ (30)

The same quasi invariance (now with a space time dependent boost parameter) is
reproduced provided the transformation laws of the newly introduced fields are
given by,

$\displaystyle\delta_{0}\Sigma_{0}{}^{0}$ $\displaystyle=$
$\displaystyle\Sigma_{0}{}^{\nu}\partial_{\nu}\zeta^{0}$
$\displaystyle\delta_{0}\Sigma_{0}{}^{k}$ $\displaystyle=$
$\displaystyle\Sigma_{0}{}^{\nu}\partial_{\nu}\zeta^{k}+u^{a}\Sigma_{a}{}^{k}$
$\displaystyle\delta_{0}\Sigma_{a}{}^{k}$ $\displaystyle=$
$\displaystyle\Sigma_{a}{}^{\nu}\partial_{\nu}\zeta^{k}-\omega_{a}{}^{b}\Sigma_{b}{}^{k}$
$\displaystyle\delta_{0}B_{\mu}$ $\displaystyle=$
$\displaystyle-\partial_{\mu}\zeta^{\nu}B_{\nu}-\frac{1}{2}\partial_{\mu}\omega^{ab}\sigma_{ab}+\partial_{\mu}u^{a}{x_{a}}$
(31)

we can easily see that the transformation relations (31) can be derived from
(29) and (10) [34] just by putting the convective term to be zero. The reason
for this is not difficult to guess. For the rela5tvistic spin less particle
and Bosonic string the parameters that describe the location of the object is
independent of choice of any coordinates. Hence when we implement spacetime
transformations the parameters do not change. The information we have about
the dynamics of the model dictate how many terms would appear in the
expression of the curved model. Thus every single new model is not going to
change the geometry. Rather it is the dynamics of the model that would
determine how many independent symmetry elements of the geometry will be
required in the model building.

##  7 Conclusion

We have discussed the formulation of the non - relativistic string coupled
with non - relativistic gravity in Newton \- Cartan manifold . This result is
in itself of paramount importance because working with a different gauging
approach it has been definitely concluded that the NC geometry breaks down to
propagate non - relativistic string. Of course we are in no position to
comment to the difference, we can definitely conclude that the gauging of
symmetry approach implemented by the Galilean gauge theory can give a string
action in the NC space time consistently.

The Galilean gauge theory is based on what may be called little figuratively —
the gauging the symmetries approach. What this means is very simple. The
target is to couple a model with Galilean symmetry in flt space time with
curved background. Since the metric structure is degenerate, coupling via the
metrics is not useful. this suggests a vielbein formulation of the theory. GGT
provides this first order formulation . Not only that, it contains the
symmetries of the model in the tangent plane inbuilt in it. No wonder, it has
been proved as a natural theory of non - relativistic diffeomorphism
invariance. The results for the string presented here completes the glorious
journey of GGT so to say [10], [11],[12], [33],[34].

## References

  * [1] C. Misner, K. Tnorn and J. Wheeler,Gravitation (Freeman, 1973)
  * [2] P. Havas, Rev. Mod. Phys. 36,(1964),938.
  * [3] J. L. Anderson: Principles of Relativity Physics, Academic Press, London 1967.
  * [4] A. Trautman, “Theories of Space, Time and Gravitation” in Lectures on General Relativity, S. Deser and K.W. Ford, eds., Prentice-hall, Englewood Cliffs, 1965.
  * [5] H.P. Künzle, Ann. Inst. Henry Poincaré, 42,(1972),337.
  * [6] K. Kuchař, Phys. Rev. ,22D,6,(1980),1285.
  * [7] G. Dautcourt: “Die Newtonske Gravitationstheorie als Strenger Grenzfall der Allgemeinen Relativitätheorie”, Acta Phys. Pol. 25,5,(1964),637.
  * [8] J. Ehlers: “On Limit Relations between, and Approximative Explanations of, Physical Theories”, in Logic, Methodology and Philosophy of Science, VII, B. Marcus et al., eds., Elsiever, Amsterdam (1986),405.
  * [9] D.B. Malament: “Newtonian Gravity, Limits, and the Geometry of Space”, in From Quarks to Quasars, ed. R.G. Colodny, Pittsburgh University Press, Pittsburgh (1986)p.181.
  * [10] R. Banerjee, A. Mitra and P. Mukherjee, Phys. Lett. B 737, 369 (2014).
  * [11] R. Banerjee, A. Mitra and P. Mukherjee, Phys. Rev. D 91, no. 8, 084021 (2015) doi:10.1103/PhysRevD.91.084021 [arXiv:1501.05468 [gr-qc]].
  * [12] R. Banerjee, A. Mitra, P. Mukherjee, Class. Quantum Grav. 32, 045010 (2015).
  * [13] D.T. Son and M. Wingate, Annals. of. Physics. 321, 197-224 (2006).
  * [14] B. Grinstein and S. Pal, Phys. Rev. D 97, no. 12, 125006 (2018)
  * [15] Eric Bergshoeff, Jaume Gomis and Ziqi Yan, arXiv:1806.06071 [hep-th].
  * [16] M. Geracie, arXiv:1611.01198 [hep-th].
  * [17] A. Jain, Phys. Rev. D 93, no. 6, 065007 (2016)
  * [18] A. Mitra, Int. J. Mod. Phys. A 32, no. 36, 1750206 (2017)
  * [19] K. Morand, arXiv:1811.12681 [hep-th].
  * [20] J. Read and N. J. Teh, Class. Quant. Grav. 35, no. 18, 18LT01 (2018)
  * [21] G. K. Karananas, doi:10.5075/epfl-thesis-7173
  * [22] E. Cartan, Sur les variétés à connexion affine et la théorie de la relativité généralisée. (première partie), Annales Sci.Ecole Norm.Sup. 40 (1923) 325–412.
  * [23] E. Cartan, Sur les variétés à connexion affine et la théorie de la relativité généralisée. (première partie) (Suite)., Annales Sci.Ecole Norm.Sup. 41 (1924) 1–25.
  * [24] S Sachdev, Quantum Phase Transition , cambridge university Press (2011) .
  * [25] K, Jensen, SciPost Phys. 5, no. 1, 011 (2018)
  * [26] R. Andringa, E. A. Bergshoeff, S. Panda and M. de Roo, Class. Quantum Grav. 28 105011 (2011).
  * [27] R. Banerjee and P. Mukherjee, Nucl. Phys. B 938, 1 (2019)
  * [28] R. Banerjee, P. Mukherjee, Phys. Rev. D93, 085020 (2016).
  * [29] R. Banerjee and P. Mukherjee, Class. Quant. Grav. 33, no. 22, 225013 (2016) doi:10.1088/0264-9381/33/22/225013 [arXiv:1604.06893 [gr-qc]].
  * [30] R. Banerjee and P. Mukherjee, Phys. Lett. B 778, 303 (2018) doi:10.1016/j.physletb.2018.01.033 [arXiv:1710.10882 [gr-qc]].
  * [31] R. Banerjee and P. Mukherjee, Nucl. Phys. B 938, 1 (2019) doi:10.1016/j.nuclphysb.2018.11.002 [arXiv:1801.08373 [gr-qc]].
  * [32] P. Mukherjee and A. Sattar, Phys. Rev. D 99, no. 8, 084038 (2019)
  * [33] R. Banerjee and P. Mukherjee, Phys. Lett. B 797, 134834 (2019) doi:10.1016/j.physletb.2019.134834 [arXiv:1907.01508 [gr-qc]].
  * [34] R. Banerjee and P. Mukherjee, Phys. Rev. D 101, no.12, 126013 (2020) doi:10.1103/PhysRevD.101.126013 [arXiv:1910.01452 [gr-qc]].
  * [35] R. Banerjee and P. Mukherjee, [arXiv:2003.05693 [gr-qc]].
  * [36] R. Banerjee , P. Mukherjee and S. Moinuddin, PhysRevD.103.046020 (2021) [arXiv:2008.00814 [gr-qc]].


93A30, 93C05, 93D09

# Parameterized Interpolation of Passive Systems

Peter Benner<EMAIL_ADDRESS>Pawan Goyal  goyalp@mpi-
magdeburg.mpg.de Paul Van Dooren<EMAIL_ADDRESS>
###### Abstract

We study the tangential interpolation problem for a passive transfer function
in standard state-space form. We derive new interpolation conditions based on
the computation of a deflating subspace associated with a selection of
spectral zeros of a parameterized para-Hermitian transfer function. We show
that this technique improves the robustness of the low order model and that it
can also be applied to non-passive systems, provided they have sufficiently
many spectral zeros in the open right half plane. We analyze the accuracy
needed for the computation of the deflating subspace, in order to still have a
passive lower order model and we derive a novel selection procedure of
spectral zeros in order to obtain low order models with a small approximation
error.

###### keywords:

Tangential interpolation, passive systems, passivity radius, robustness

  * •

We investigate an tangential interpolation problem for passive systems.

  * •

We propose interpolation conditions based on a deflating subspace.

  * •

We discuss a construction of passive low-order models using the deflating
subspace associated with a selection of spectral zeros.

##  1 Introduction

We consider linear and finite dimensional dynamical systems that are
_passive_. We restrict ourselves to continuous-time systems that can be
represented in standard state-space form with real coefficients and real
inputs, outputs and states :

$\begin{array}[]{rcl}\dot{x}(t)&=&Ax(t)+Bu(t),\ x(0)=0,\\\
y(t)&=&Cx(t)+Du(t).\end{array}$ (1.1)

Denoting real and complex $n$-vectors ($n\times m$ matrices) by
$\mathbb{R}^{n}$, $\mathbb{C}^{n}$ ($\mathbb{R}^{n\times m}$,
$\mathbb{C}^{n\times m}$), respectively, then $u:\mathbb{R}\to\mathbb{R}^{m}$,
$x:\mathbb{R}\to\mathbb{R}^{n}$, and $y:\mathbb{R}\to\mathbb{R}^{m}$ are
vector-valued functions denoting the _input_ , _state_ , and _output_ of the
system, and the coefficient matrices satisfy $A\in\mathbb{R}^{n\times n}$,
$B\in\mathbb{R}^{n\times m}$, $C\in\mathbb{R}^{m\times n}$, and
$D\in\mathbb{R}^{m\times m}$.

Model reduction of such systems has been a major research topic for the last
three decades and led to a wealth of different approaches, as illustrated in
several survey volumes [2], [5], [6]. One of the proposed approaches is based
on tangential interpolation [10], [3]. This technique was originally developed
for arbitrary types of rational transfer functions [10], but an important
drawback is that some critical properties – such as stability or passivity –
are not easy to satisfy and require a careful selection of interpolation
conditions. It was shown in [3] that when using spectral zeros of a given
transfer function as interpolation conditions, then one does preserve
passivity in the reduced-order model, at least for the single-input/single-
output case. This was extended by Sorensen in [14] to the multi-input/multi-
output case by making use of deflating subspace calculations and the Kalman-
Yakubovich–Popov conditions for passivity. Numerical and structure-preserving
algorithms to compute reduced-order models based on this approach are
suggested in [8]. The link between both methods was later on pointed out by
Fanizza et al. [9], who also make the connection to the so-called covariance
extension problem. In the present paper, we further extend the approach of
Sorensen by applying it to a class of systems that are parameterized by a
scalar parameter. The new contributions of this paper are threefold :

  1. 1.

we show that we can apply the deflating subspace idea to a class of
parameterized systems, which improves the robustness of the reduced-order
system by increasing its passivity radius,

  2. 2.

we derive a novel selection technique of the subset of spectral zeros used for
model reduction, which attempts to minimize the approximation error, and

  3. 3.

we show that the method can be applied to non-passive systems and still
constructs passive lower order models, under certain conditions.

Because of the last property, we give in this paper a new derivation of
Sorensen’s results in order to show that it may apply also to non-passive
systems.

The notation used in the paper is as follows. The Hermitian (or conjugate)
transpose (transpose) of a vector or matrix $V$ is denoted by $V^{\mathsf{H}}$
($V^{\mathsf{T}}$) and the identity matrix is denoted by $I_{n}$ or $I$ if the
dimension is clear. We require that input and output dimensions are equal to
$m$ since we want to interpolate with passive transfer functions. Throughout
this article we will use the following notation. We denote the set of
symmetric matrices in $\mathbb{R}^{n\times n}$ by ${\mathbb{S}_{n}}$. Positive
definiteness (semi-definiteness) of $M\in{\mathbb{S}_{n}}$ is denoted by
$M\succ 0$ ($M\succeq 0$). In Section 2, we recall the properties of passive
and of port-Hamiltonian systems in order to define the robustness measure
known as the passivity radius. In Section 3 we then recall the results of
Sorensen on interpolation in spectral zeros via deflating subspace
calculations. This is then extended in Section 4 to the novel technique of
parameterized passive interpolation. The selection technique to find
appropriate spectral zeros to minimize the approximation error is described in
Section 5, and the resulting robustness property is briefly described in
Section 6. We give numerical experiments that illustrate our new method and
its properties in Section 7 and give come concluding remarks in Section 8.

##  2 Passive systems and port-Hamiltonian realizations

Passive systems are well studied in the continuous-time case. We briefly
recall some important properties following [16], and refer to the literature
for a more detailed survey. We consider continuous-time systems with a
rational transfer matrix $Z(s)$ and define the following rational matrix
function of $s\in\mathbb{C}$ :

$\Phi(s):=Z^{\mathsf{T}}(-s)+Z(s),$

which is said to be _para-Hermitian_ since $\Phi(-s)^{\mathsf{T}}=\Phi(s)$. It
therefore coincides with two times the Hermitian part of $Z(s)$ on the
$\imath\omega$ axis:

$\Phi(\imath\omega)=[Z(\imath\omega)]{{}^{\mathsf{H}}}+Z(\imath\omega).$

######  Definition 1.

The rational transfer function $Z(s)$ is called _strictly positive-real_ if
$\Phi(\imath\omega)\succ 0$ for all $\omega\in\mathbb{R}$ and it is called
_positive-real_ if $\Phi(\imath\omega)\succeq 0$ for all $\
\omega\in\mathbb{R}$.

The transfer function $Z(s)$ is called _asymptotically stable_ if the
eigenvalues of $A$ are in the open left half plane, and it is called _stable_
if the eigenvalues of $A$ are in the closed left half plane, with any
eigenvalues occurring on the imaginary axis being semi-simple.

The transfer function $Z(s)$ is called _strictly passive_ if it is strictly
positive-real and asymptotically stable and it is called _passive_ if it is
positive real and stable with polar residues that are Hermitian and positive
semi-definite for every pole on the imaginary axis.

######  Remark 1.

In the classical circuit theory literature the notion of positive realness is
phrased differently and it implies stability. For rational transfer functions
it is therefore equivalent to passivity [1]. In this paper, though, we will
use the above modified definition of positive realness.

In this paper we focus on systems that are strictly passive, which implies
that the transfer matrix has no infinite or imaginary axis poles, and hence is
proper. Moreover, $\Phi(\imath\omega)\succ 0$ at $\omega=\infty$ implies that
$D^{\mathsf{T}}+D\succ 0$ and that $\Phi(s)$ is regular. We will see that this
restriction simplifies our discussion significantly. This is also a reasonable
restriction because passive systems can be viewed as limiting cases of
strictly passive systems.

Since the transfer function is proper, we can represent it in standard state-
space form $Z(s)=C(sI_{n}-A)^{-1}B+D$ and we will assume throughout the paper
that this realization is minimal (i.e. controllable and observable). We can
associate with $\Phi(s)$ a system matrix $S(s)$ which is a generalized state-
space realization of $\Phi(s)$ :

${S(s):=\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\\
A^{\mathsf{T}}+sI_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right].}$ (2.1)

If the quadruple ${\mathcal{M}}:=\left\\{A,B,C,D\right\\}$ is a minimal
realization of a strictly passive transfer function $Z(s)$ of McMillan degree
$n$, then $S(s)$ is a minimal realization (in generalized state-space form) of
$\Phi(s)$. This transfer function has indeed degree $2n$ since $Z(s)$ and
$Z^{\mathsf{T}}(-s)$ have no common poles because of the assumption that
$Z(s)$ is asymptotically stable. Since $D^{\mathsf{T}}+D$ is nonsingular, the
$2n$ finite eigenvalues of the pencil $S(s)$ are then the so-called spectral
zeros of the strictly passive transfer function $Z(s)$.

We can apply the following congruence transformation to $S(s)$, using a
symmetric matrix $X$ :

$S_{X}(s)=\left[\begin{array}[]{cc|c}I_{n}&0&0\\\ -X&I_{n}&0\\\ \hline\cr
0&0&I_{m}\end{array}\right]\\!S(s)\\!\left[\begin{array}[]{cc|c}I_{n}&-X&0\\\
0&I_{n}&0\\\ \hline\cr 0&0&I_{m}\end{array}\right]$
${=\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\\
A^{\mathsf{T}}+sI_{n}&-A^{\mathsf{T}}X-XA&C^{\mathsf{T}}-XB\\\ \hline\cr
B^{\mathsf{T}}&C-B^{\mathsf{T}}X&D^{\mathsf{T}}+D\end{array}\right]}$ (2.2)

without affecting the transfer function $\Phi(s)$ of this system matrix. If
the following submatrix of $S_{X}(s)$

${\mathcal{W}}(X,{\mathcal{M}}):=\left[\begin{array}[]{cc}-A^{\mathsf{T}}X-X\,A&C^{\mathsf{T}}-X\,B\\\
C-B^{\mathsf{T}}X&D^{\mathsf{T}}+D\end{array}\right]$ (2.3)

is positive semi-definite, then it can be factored as indicated below

${{\mathcal{W}}(X,{\mathcal{M}})=\left[\begin{array}[]{c}C_{G}^{\mathsf{T}}\\\
D_{G}^{\mathsf{T}}\end{array}\right]\left[\begin{array}[]{cc}C_{G}&D_{G}\end{array}\right],}$

from which it follows that

${S_{X}(s)=\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\\
A^{\mathsf{T}}+sI_{n}&C_{G}^{\mathsf{T}}C_{G}&C_{G}^{\mathsf{T}}D_{G}\\\
\hline\cr
B^{\mathsf{T}}&D_{G}^{\mathsf{T}}C_{G}&D_{G}^{\mathsf{T}}D_{G}\end{array}\right]}$

and hence $G(s):=C_{G}(s\,I_{n}-A)^{-1}B+D_{G}$ is the right factor of the
symmetric factorization $\Phi(s)=G^{\mathsf{T}}(-s)G(s)$. This then implies
that $\Phi(s)$ is positive semi-definite on the $\imath\omega$ axis. Moreover,
if $A$ is assumed to be asymptotically stable, then the matrix $X$ in (2.3)
must be positive definite. This discussion is an intuitive explanation of the
following result, a proof of which can be found in [16], [11].

######  Theorem 1.

Let ${\mathcal{M}}:=\left\\{A,B,C,D\right\\}$ be a minimal realization of a
proper rational transfer function $Z(s)$ and let
${\mathcal{W}}(X,{\mathcal{M}})$ be the associated matrix defined in (2.3).
Then $Z(s)$
(i) is positive real if and only if there exists a real symmetric matrix
$X\in{\mathcal{S}}$ such that

${\mathcal{W}}(X,{\mathcal{M}})\succeq 0,$ (2.4)

(ii) is passive if and only if there exists a real symmetric matrix
$X\in{\mathcal{S}}$ such that

${\mathcal{W}}(X,{\mathcal{M}})\succeq 0,\quad X\succ 0,$ (2.5)

and (iii) is strictly passive if and only if there exists a real symmetric
matrix $X\in{\mathcal{S}}$ such that

${\mathcal{W}}(X,{\mathcal{M}})\succ 0,\quad X\succ 0.$ (2.6)

The Linear Matrix Inequality (LMI) given is (2.5) is also known as the Kalman-
Yakubovich-Popov condition for passivity. In the sequel we will make use of
the solution sets of these inequalities :

$\displaystyle{\mathbb{X}^{\raisebox{1.39998pt}{{$>$}}}}:=\left\\{X\in{\mathcal{S}}\left|{\mathcal{W}}(X,{\mathcal{M}})\succeq
0,\ X\succ 0\right.\right\\},$ (2.7a)
$\displaystyle{\mathbb{X}^{\raisebox{1.39998pt}{{$\gg$}}}}:=\left\\{X\in{\mathcal{S}}\left|{\mathcal{W}}(X,{\mathcal{M}})\succ
0,\ X\succ 0\right.\right\\}.$ (2.7b)

######  Definition 2.

Every solution of the LMI (2.7a) is called a certificate for the passivity of
the model ${\mathcal{M}}$ and every solution of the LMI (2.7b) is called a
certificate for the strict passivity of the model ${\mathcal{M}}$.

If $D^{\mathsf{T}}+D$ is invertible, then the solutions in
${\mathbb{X}^{\raisebox{1.39998pt}{{$>$}}}}$ where
${\mathcal{W}}(X,{\mathcal{M}})$ is of minimum rank, are those for which
$\operatorname{rank}{\mathcal{W}}(X,{\mathcal{M}})=\operatorname{rank}(D^{\mathsf{T}}+D)=m$,
which is the case if and only if the Schur complement of $D^{\mathsf{T}}+D$ in
${\mathcal{W}}(X,{\mathcal{M}})$ is zero. This Schur complement is associated
with the continuous-time _algebraic Riccati equation (ARE)_

$\mathsf{Ricc}(X):=-XA-A^{\mathsf{T}}X-(C^{\mathsf{T}}-XB)(D^{\mathsf{T}}+D)^{-1}(C-B^{\mathsf{T}}X)=0.$
(2.8)

Each symmetric solution $X$ to (2.8) yields a spectral factorization
$\Phi(s)=G^{\mathsf{T}}(-s)G(s)$ where $G(s)$ is $m\times m$ and regular.
Therefore, the spectral zeros of $\Phi(s)$ are the union of the zeros of
$G(s)$ and of $G^{\mathsf{T}}(-s)$. The matrix $X$ also corresponds to an
invariant subspace spanned by the columns of
$U:=\left[\begin{array}[]{cc}I_{n}\\\ -X\end{array}\right]$ that remains
invariant under multiplication with the _Hamiltonian matrix_

$H:=\left[\begin{array}[]{cc}A-B(D^{\mathsf{T}}+D)^{-1}C&-B(D^{\mathsf{T}}+D)^{-1}B^{\mathsf{T}}\\\
C^{\mathsf{T}}(D^{\mathsf{T}}+D)^{-1}C&-(A-B(D^{\mathsf{T}}+D)^{-1}C)^{\mathsf{T}}\end{array}\right],$
(2.9)

i. e. $U$ satisfies $HU=UA_{F}$ for a matrix $A_{F}=A-BF$ with
$F:=(D^{\mathsf{T}}+D)^{-1}(C-B^{\mathsf{T}}X)$. We point out here that the
solutions $X$ of the Riccati equations are certificates for the passivity of
the model ${\mathcal{M}}$, but not for its strict passivity. We will see that
this distinction plays an important role in the sequel. It is also shown in
[16] that for a minimal model ${\mathcal{M}}$, the set of solutions $X$ of the
Riccati equation (2.8) has two extremal solutions $X_{-}$ and $X_{+}$ such
that all other certificates $X\in{\mathbb{X}^{\raisebox{1.39998pt}{{$>$}}}}$
satisfy $X_{-}\preceq X\preceq X_{+}$.

We now give a brief introduction to special realizations of passive systems,
known as port-Hamiltonian system models.

######  Definition 3.

A linear time-invariant _port-Hamiltonian (pH) system model_ of a proper
transfer function, has the standard state-space form

$\begin{array}[]{rcl}\dot{x}&=&(J-R)Qx+(G-P)u,\\\
y&=&(G+P)^{\mathsf{T}}Qx+(N+S)u,\end{array}$ (2.10)

where the system matrices satisfy the symmetry conditions

${\mathcal{V}}:=\left[\begin{array}[]{cccc}-J&-G\\\
G^{\mathsf{T}}&N\end{array}\right]=-\mathcal{V}^{\mathsf{T}},\
{\mathcal{W}}:=\left[\begin{array}[]{cccc}R&P\\\
P^{\mathsf{T}}&S\end{array}\right]=\mathcal{W}^{\mathsf{T}}\succeq 0,\
Q=Q^{\mathsf{T}}\succeq 0.$ (2.11)

Port-Hamiltonian systems were introduced from a different point of view [15],
but they are also known to be passive. If the model is strictly passive then
$A$ and hence also $Q$ are both invertible. We can then choose $X=Q$ as
certificate to show that the model
${\mathcal{M}}:=\left\\{(J-R)Q,G-P,(G+P)^{\mathsf{T}}Q,N+S\right\\}$ satisfies
the KYP condition. Conversely, let ${\mathcal{M}}:=\left\\{A,B,C,D\right\\}$
be a state-space model satisfying the KYP condition (2.5) with a given $X\succ
0$. Then it can always be put in port-Hamiltonian form, as indicated in [4].
We can use a symmetric factorization $X=T^{\mathsf{T}}T$, which implies the
invertibility of $T$, and define a new realization

$\\{A_{T},B_{T},C_{T},D_{T}\\}:=\\{TAT^{-1},TB,CT^{-1},D\\}$

so that

$\displaystyle\left[\begin{array}[]{cccc}T^{-\mathsf{T}}&0\\\
0&I_{m}\end{array}\right]\left[\begin{array}[]{cccc}-A^{\mathsf{T}}X-XA&C^{\mathsf{T}}-XB\\\
C-B^{\mathsf{T}}X&D^{\mathsf{T}}+D\end{array}\right]\left[\begin{array}[]{cccc}T^{-1}&0\\\
0&I_{m}\end{array}\right]$ (2.18)
$\displaystyle\qquad=\left[\begin{array}[]{cccc}-A_{T}&-B_{T}\\\
C_{T}&D_{T}\end{array}\right]+\left[\begin{array}[]{cccc}-A_{T}^{\mathsf{T}}&C^{\mathsf{T}}_{T}\\\
-B^{\mathsf{T}}_{T}&D^{\mathsf{T}}_{T}\end{array}\right]\succeq 0.$ (2.23)

We can then use the symmetric and skew-symmetric part of the matrix

${\mathcal{S}}:=\left[\begin{array}[]{cccc}-A_{T}&-B_{T}\\\
C_{T}&D_{T}\end{array}\right]$

to define the coefficients of a pH representation via

${\mathcal{V}}:=\left[\begin{array}[]{cccc}J&G\\\
-G^{\mathsf{T}}&N\end{array}\right]:=\frac{{{\mathcal{S}}}-{{\mathcal{S}}}^{\mathsf{T}}}{2},\quad{\mathcal{W}}:=\left[\begin{array}[]{cccc}R&P\\\
P^{\mathsf{T}}&S\end{array}\right]:=\frac{{{\mathcal{S}}}+{{\mathcal{S}}}^{\mathsf{T}}}{2}\succeq
0.$

This construction yields $Q=I_{n}$ because of the chosen factorization
$X=T^{\mathsf{T}}T$. A system with such a representation will be called a
normalized port-Hamiltonian system. This shows that proper port-Hamiltonian
systems are nothing but passive systems described in an appropriate coordinate
system. On the other hand, the passivity radius of a normalized port-
Hamiltonian system has good robustness properties in terms of its so-called
passivity radius, as is shown below.

######  Definition 4.

The _passivity radius_ $\rho_{\mathcal{M}}$ of a passive model
${\mathcal{M}}:=\\{A,B,C,D\\}$ is the smallest perturbation
$\Delta_{\mathcal{M}}:=\\{\Delta_{A},\Delta_{B},\Delta_{C},\Delta_{D}\\}$
which can make the model ${\mathcal{M}}+\Delta_{\mathcal{M}}$ loose its
passivity.

Therefore, if the perturbation $\Delta_{\mathcal{M}}$ is measured by

$\|\Delta_{\mathcal{M}}\|:=\left\|\left[\begin{array}[]{cc}\Delta_{A}&\Delta_{B}\\\
\Delta_{C}&\Delta_{D}\end{array}\right]\right\|_{2},\quad\mathrm{or}\quad\left\|\left[\begin{array}[]{cc}\Delta_{A}&\Delta_{B}\\\
\Delta_{C}&\Delta_{D}\end{array}\right]\right\|_{F},$

then, for a given certificate
$X\in{\mathbb{X}^{\raisebox{1.39998pt}{{$\gg$}}}}$, the passivity condition
${\mathcal{W}}(X,{\mathcal{M}}+\Delta_{\mathcal{M}})\succeq 0$ for all
perturbed systems ${\mathcal{M}}+\Delta_{\mathcal{M}}$ becomes just a linear
matrix inequality in $\Delta_{\mathcal{M}}$, and hence yields a computable
lower bound for $\rho_{\mathcal{M}}$, which is called the $X$-passivity radius
$\rho_{\mathcal{M}}(X)$ :

$\rho_{\mathcal{M}}(X):=\inf_{\Delta_{\mathcal{M}}\in\mathbb{C}^{n+m,n+m}}\left\\{\|\Delta_{\mathcal{M}}\|\;|\;\det{\mathcal{W}}(X,{\mathcal{M}}+\Delta_{\mathcal{M}})=0\right\\}\leq\rho_{\mathcal{M}}.$

It follows (see [13]) that $\rho_{\mathcal{M}}$ is the supremum of these lower
bounds over all certificates
$X\in{\mathbb{X}^{\raisebox{1.39998pt}{{$\gg$}}}}$ :

$\rho_{{\mathcal{M}}}:=\sup_{X\in{\mathbb{X}^{\raisebox{0.99998pt}{{$\gg$}}}}}\rho_{{\mathcal{M}}}(X).$
(2.24)

The following theorem, proven in [13], shows that normalized port-Hamiltonian
systems have an $X$-passivity radius that is at least as good as the
corresponding non-normalized system.

######  Theorem 2.

Let ${\mathcal{M}}=\\{A,B,C,D\\}$ be a model of a strictly passive transfer
function $Z(s)$. Then for every certificate
$X\in{\mathbb{X}^{\raisebox{1.39998pt}{{$\gg$}}}}$, we can construct a
normalized port-Hamiltonian system

${\mathcal{M}}_{T}:=\\{A_{T},B_{T},C_{T},D_{T}\\}=\left\\{J-R,G-P,(G+P)^{\mathsf{T}},N+S\right\\}$

using a factorization $X=T^{\mathsf{T}}T$. The $X$-passivity radii
$\rho_{\mathcal{M}}(X)$ and $\rho_{{\mathcal{M}}_{T}}(I)$ of these two models
satisfy

$\rho_{\mathcal{M}}(X)\leq\rho_{{\mathcal{M}}_{T}}(I)=\lambda_{\min}({\mathcal{W}}).$

The optimal passivity radius for all possible models for $Z(s)$ must therefore
be attained by a normalized port-Hamiltonian model. The following theorem
indicates that there is such a normalized port-Hamiltonian system with optimal
passivity radius and that it corresponds to a certificate $X$ for a family of
passive systems, parameterized by the real parameter $\xi$ :

${\mathcal{M}}_{\xi}:=\\{A+\frac{\xi}{2}I_{n},B,C,D-\frac{\xi}{2}I_{m}\\},\quad
Z_{\xi}(s):=C((s-\xi/2)I_{n}-A)^{-1}B+(D-\xi I_{m}/2).$ (2.25)

######  Theorem 3.

Let $Z(s)$ be a given strictly passive transfer function, then there exists a
port-Hamiltonian system model ${\mathcal{M}}$ of $Z(s)$ with the largest
possible passivity radius, and it corresponds to a common certificate $X$ for
all the transfer functions $Z_{\xi}(s)$ that are strictly passive where
$0<\xi<\Xi$ and $\Xi$ is the smallest positive number such that $Z_{\Xi}(s)$
is passive, but not strictly passive.

It was shown in [13] that the calculation of $\Xi$ is a two dimensional
optimization problem that can be solved efficiently. Once the value of $\Xi$
is known, one can find a certificate $X\succ 0$ for the LMI
${\mathcal{W}}(X,{\mathcal{M}}_{\Xi})\succeq 0$ using a Riccati equation
approach or the corresponding generalized eigenvalue problem

${S_{\Xi}(s):=\left[\begin{array}[]{cc|c}0&A+\Xi I_{n}/2-sI_{n}&B\\\
A^{\mathsf{T}}+\Xi I_{n}/2+sI_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D-\Xi I_{m}\end{array}\right].}$

That certificate is then valid for the family of LMIs
${\mathcal{W}}(X,{\mathcal{M}}_{\xi})\succ 0$ for $0<\xi<\Xi$, and indicating
that the transfer functions $Z_{\xi}(s):=C((s-\xi/2)I_{n}-A)^{-1}B+(D-\xi
I_{m}/2)$ are all strictly passive. We will use this fact later on to propose
a family of interpolation conditions of low order passive transfer functions
approximating a high order one.

##  3 Passive interpolation using spectral zeros

We rederive here the technique developed by Sorensen [14] and Antoulas [3] for
the construction of a degree $\hat{n}$ passive system $\hat{Z}(s)$
approximating a given passive system $Z(s)$ of McMillan degree $n\geq\hat{n}$,
via interpolation in a set of so-called spectral zeros. But we relax the
conditions imposed on the transfer function $Z(s)$, because we will need this
in the next section. Our proof is based on Sorenson’s construction, but it
applies also to non-passive systems $Z(s)$.

######  Theorem 4.

Let ${\mathcal{M}}:=\\{A,B,C,D\\}$ be a minimal model of an $m\times m$
transfer function $Z(s)$ and let $S(s)$ be the system matrix of $\Phi(s)$.
Assume that $D+D^{\mathsf{T}}\succ 0$ and that we are then given a basis for
an $\hat{n}$ dimensional deflating subspace of $S(s)$ satisfying

$\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\\
A^{\mathsf{T}}+sI_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right]\left[\begin{array}[]{c}U\\\
V\\\ \hline\cr W\end{array}\right]=\left[\begin{array}[]{c}V\\\ -U\\\
\hline\cr 0\end{array}\right](R-sI_{\hat{n}}),$ (3.1)

where the spectrum of $R$ lies in the open right half plane. Then
$\hat{X}:=-U^{\mathsf{T}}V$ is symmetric. If, moreover, $\hat{X}$ is
invertible, then the reduced-order transfer function $\hat{Z}(s)$ of the
projected system model

$\hat{\mathcal{M}}:=\\{\hat{A},\hat{B},\hat{C},\hat{D}\\}=\\{(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}AV,(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}B,CV,D\\}$

is strictly positive real, and if $\hat{X}$ is also positive definite, then
$\hat{Z}(s)$ is strictly passive.

###### Proof.

The proof uses several arguments given in [14] for the more restrictive
problem of a strictly passive transfer function $Z(s)$. The symmetry of
$\hat{X}:=-U^{\mathsf{T}}V$ follows from the following equation, obtained from
multiplying (3.1) on the left with
$\left[U^{\mathsf{T}}\;V^{\mathsf{T}}\;W^{\mathsf{T}}\right]$ :

$\left[\begin{array}[]{cc|c}U^{\mathsf{T}}\\!&\\!V^{\mathsf{T}}\\!&\\!W^{\mathsf{T}}\\!\end{array}\right]\\!\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\\
A^{\mathsf{T}}+sI_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right]\left[\begin{array}[]{c}U\\\
V\\\ \hline\cr
W\end{array}\right]=(U^{\mathsf{T}}V-V^{\mathsf{T}}U)(R-sI_{\hat{n}}).$

Since the left hand side is para-Hermitian, the right hand side must also be
para-Hermitian, which implies that

$(U^{\mathsf{T}}V-V^{\mathsf{T}}U)(R-sI_{\hat{n}})=(R^{\mathsf{T}}+sI_{\hat{n}})(V^{\mathsf{T}}U-U^{\mathsf{T}}V),$

and finally,

$(U^{\mathsf{T}}V-V^{\mathsf{T}}U)R+R^{\mathsf{T}}(U^{\mathsf{T}}V-V^{\mathsf{T}}U)=0.$

Since $R$ has all its eigenvalues in the right half plane, the matrix
$(U^{\mathsf{T}}V-V^{\mathsf{T}}U)$ must be zero, which implies that
$\hat{X}:=-U^{\mathsf{T}}V$ is symmetric. Using the different rows of (3.1)
one obtains

$\displaystyle AV+BW$ $\displaystyle=$ $\displaystyle VR$ (3.2) $\displaystyle
A^{\mathsf{T}}U+C^{\mathsf{T}}W$ $\displaystyle=$ $\displaystyle-UR$ (3.3)
$\displaystyle B^{\mathsf{T}}U+CV+(D^{\mathsf{T}}+D)W$ $\displaystyle=$
$\displaystyle 0$ (3.4)

and from equations (3.2), (3.3) and the symmetry of $\hat{X}$, it follows that

$(U^{\mathsf{T}}AV+V^{\mathsf{T}}A^{\mathsf{T}}U)+(U^{\mathsf{T}}B+V^{\mathsf{T}}C)W=0.$
(3.5)

If $\hat{X}=-U^{\mathsf{T}}V$ is invertible, we can construct the reduced-
order system model

$\hat{\mathcal{M}}:=\\{\hat{A},\hat{B},\hat{C},\hat{D}\\}=\\{(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}AV,(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}B,CV,D\\}$
(3.6)

with transfer function $\hat{Z}(s)$, and it then follows from (3.4) and (3.5)
that

${\mathcal{W}}(\hat{X},\hat{\mathcal{M}}):=\left[\begin{array}[]{cc}-\hat{A}^{\mathsf{T}}\hat{X}-\hat{X}\,\hat{A}&\hat{C}^{\mathsf{T}}-\hat{X}\,\hat{B}\\\
\hat{C}-\hat{B}^{\mathsf{T}}\hat{X}&\hat{D}^{\mathsf{T}}+\hat{D}\end{array}\right]=\left[\begin{array}[]{c}\\!-W^{\mathsf{T}}\\!\\\
I_{m}\end{array}\right](\hat{D}^{\mathsf{T}}+\hat{D})\left[\begin{array}[]{cc}-W&I_{m}\\!\end{array}\right]\succeq
0,$ (3.7)

since
${\mathcal{W}}(\hat{X},\hat{\mathcal{M}})\left[\begin{array}[]{c}I_{\hat{n}}\\\
W\end{array}\right]=0$ and $\hat{D}^{\mathsf{T}}+\hat{D}\succ 0$. It then
follows from Theorem 1 that the function
$\hat{\Phi}(s)=\hat{Z}^{\mathsf{T}}(-s)+\hat{Z}(s)$ is non-negative on the
imaginary axis, and hence that $\hat{Z}(s)$ is positive real. Moreover, (3.7)
implies that

$\hat{\Phi}(s)=\hat{G}^{\mathsf{T}}(-s)(D^{\mathsf{T}}+D)\hat{G}(s),\quad\mathrm{where}\quad
G(s)=I_{m}-W(sI_{\hat{n}}-\hat{A})^{-1}\hat{B}$

and has as zeros the eigenvalues of $R$ since (3.2) implies that
$\hat{A}+\hat{B}W=R$. Therefore, $\hat{\Phi}(s)$ has no zeros on the imaginary
axis and hence must be strictly positive real. Finally, if $\hat{X}$ is
positive definite, then $\hat{A}$ is also asymptotically stable, which means
that $\hat{Z}(s)$ is strictly passive. ∎

######  Remark 2.

In [14], Sorensen proves that if $Z(s)$ is strictly passive, then $\Phi(s)$
has $n$ spectral zeros in the right half plane and the conditional assumptions
of the above theorem always hold true for every choice of $\hat{n}\leq n$
right plane spectral zeros, implying that $\hat{X}$ is positive definite and
$\hat{Z}(s)$ is strictly passive.

######  Corollary 3.1.

The above equation (3.7) indicates that
${\mathcal{W}}(\hat{X},\hat{\mathcal{M}})$ has minimum rank, which implies
that its Schur complement is zero, and hence that $\hat{X}$ solves the Riccati
equation

$-\hat{X}\hat{A}-\hat{A}^{\mathsf{T}}\hat{X}-(\hat{C}^{\mathsf{T}}-\hat{X}\hat{B})(\hat{D}^{\mathsf{T}}+\hat{D})^{-1}(\hat{C}-\hat{B}^{\mathsf{T}}X)=0.$

Moreover, the corresponding feedback matrix
$F:=(\hat{D}^{\mathsf{T}}+\hat{D})^{-1}(\hat{C}-\hat{B}^{\mathsf{T}}X)$ is a
stabilizing feedback since ${\mathcal{W}}(\hat{X},\hat{\mathcal{M}})\succeq
0$.

It was shown in [9] that when $Z(s)$ is strictly passive and $R$ has distinct
eigenvalues (which is the generic case), then the lower order model
$\hat{\mathcal{M}}:=\\{\hat{A},\hat{B},\hat{C},\hat{D}\\}$ constructed as in
(3.6), satisfies the following tangential interpolation conditions :

$Z(\lambda_{j})Wr_{j}=\hat{Z}(\lambda_{j})Wr_{j},\;\;r^{\mathsf{T}}_{j}W^{\mathsf{T}}Z(-\lambda_{j})=r^{\mathsf{T}}_{j}W^{\mathsf{T}}\hat{Z}(-\lambda_{j}),\;\;j=1,...,\hat{n},\;\;Z(\infty)=\hat{Z}(\infty),$

where $(\lambda_{j},r_{j}),\;j=1,\ldots,\hat{n}$, is a set of self-conjugate
(eigenvalue, eigenvector) pairs of $R$. When $R$ has distinct eigenvalues,
this relates the method of Sorenson to the spectral zero interpolation
approach of Antoulas [3, 2]. If some of the eigenvalues are repeated, the
conditions imply also that derivatives at these points should match (see [10,
2]).

We give below a more complete (and simpler) proof of this connection, for the
case where $Z(s)$ satisfies the relaxed conditions of Theorem 4.

######  Theorem 5.

Let ${\mathcal{M}}:=\\{A,B,C,D\\}$ be a minimal model of an $m\times m$
transfer function $Z(s)$ and let $S(s)$ be the system matrix of $\Phi(s)$.
Assume that $D+D^{\mathsf{T}}\succ 0$ and that we are then given a basis for
an $\hat{n}$ dimensional deflating subspace of $S(s)$ satisfying

$\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\\
A^{\mathsf{T}}+sI_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right]\left[\begin{array}[]{c}U\\\
V\\\ \hline\cr W\end{array}\right]=\left[\begin{array}[]{c}V\\\ -U\\\
\hline\cr 0\end{array}\right](R-sI_{\hat{n}}),$ (3.8)

where the spectrum of $R$ lies in the open right half plane and the matrix
$\hat{X}:=-U^{\mathsf{T}}V$ is positive definite. Then the reduced-order
transfer function $\hat{Z}(s)$ of the projected system

$\hat{\mathcal{M}}:=\\{\hat{A},\hat{B},\hat{C},\hat{D}\\}=\\{(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}AV,(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}B,CV,D\\}$

is strictly positive real, and it satisfies the following tangential
interpolation conditions that define $\hat{Z}(s)$ completely :

$Z(\lambda_{j})Wr_{j}=\hat{Z}(\lambda_{j})Wr_{j},\;\;r^{\mathsf{T}}_{j}W^{\mathsf{T}}Z(-\lambda_{j})=r^{\mathsf{T}}_{j}W^{\mathsf{T}}\hat{Z}(-\lambda_{j}),\;j=1,...,\hat{n},\;\;Z(\infty)=\hat{Z}(\infty),$
(3.9)

where $(\lambda_{j},r_{j}),\;j=1,\ldots,\hat{n}$, is a set of self-conjugate
(eigenvalue, eigenvector) pairs of $R$.

###### Proof.

When multiplying the columns of (3.1) with $r_{j}$, and evaluating this at
$\lambda_{j}$, we obtain

$\left[\begin{array}[]{cc|c}0&A-\lambda_{j}I_{n}&B\\\
A^{\mathsf{T}}+\lambda_{j}I_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right]\left[\begin{array}[]{c}Ur_{j}\\\
Vr_{j}\\\ \hline\cr Wr_{j}\end{array}\right]=0,$

which implies that $Wr_{j}$ is in the kernel of the Schur complement of the
system matrix on the left :

$\Phi(\lambda_{j})Wr_{j}=\left(Z^{\mathsf{T}}(-\lambda_{j})+Z(\lambda_{j})\right)Wr_{j}=0.$
(3.10)

It follows also from (3.1) that the projected system matrix

$\left[\begin{array}[]{cc|c}0&\hat{A}-sI_{n}&\hat{B}\\\
\hat{A}^{\mathsf{T}}+sI_{n}&0&\hat{C}^{\mathsf{T}}\\\
\hline\cr\hat{B}^{\mathsf{T}}&\hat{C}&\hat{D}^{\mathsf{T}}+\hat{D}\end{array}\right]:=$
$\left[\begin{array}[]{cc|c}\\!(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}\\!&&\\\
&V^{\mathsf{T}}&\\\
\hline\cr&&I_{m}\end{array}\right]\\!\left[\begin{array}[]{cc|c}0&\\!A-sI_{n}\\!&B\\\
\\!A^{\mathsf{T}}+sI_{n}\\!&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right]\\!\left[\begin{array}[]{cc|c}\\!U(V^{\mathsf{T}}U)^{-1}\\!&&\\\
&V&\\\ \hline\cr&&I_{m}\end{array}\right]$

has $\hat{\Phi}(s):=\hat{Z}^{\mathsf{T}}(-s)+\hat{Z}(s)$ as Schur complement.
Since we have

$\left[\begin{array}[]{cc|c}U(V^{\mathsf{T}}U)^{-1}&&\\\ &V&\\\
\hline\cr&&I_{m}\end{array}\right]\left[\begin{array}[]{c}(V^{\mathsf{T}}U)r_{j}\\\
r_{j}\\\ \hline\cr Wr_{j}\end{array}\right]=\left[\begin{array}[]{c}Ur_{j}\\\
Vr_{j}\\\ \hline\cr Wr_{j}\end{array}\right]$

it follows that

$\left[\begin{array}[]{cc|c}0&\hat{A}-\lambda_{j}I_{n}&\hat{B}\\\
\hat{A}^{\mathsf{T}}+\lambda_{j}I_{n}&0&\hat{C}^{\mathsf{T}}\\\
\hline\cr\hat{B}^{\mathsf{T}}&\hat{C}&\hat{D}^{\mathsf{T}}+\hat{D}\end{array}\right]\left[\begin{array}[]{c}(V^{\mathsf{T}}U)r_{j}\\\
r_{j}\\\ \hline\cr Wr_{j}\end{array}\right]=0$

which then in turn implies that

$\hat{\Phi}(\lambda_{j})Wr_{j}=\left(\hat{Z}^{\mathsf{T}}(-\lambda_{j})+\hat{Z}(\lambda_{j})\right)Wr_{j}=0.$
(3.11)

This shows that the spectral zeros $\lambda_{j}$ and corresponding zero
directions $Wr_{j},j=1,\ldots,\hat{n}$, of $\hat{\Phi}(s)$ are a subset of
those of the original system $\Phi(s)$. To show that this also implies (3.9)
we use the same reasoning as above to obtain the equations

$\left[\begin{array}[]{cc}A-\lambda_{j}I_{n}&B\\\
C&D\end{array}\right]\left[\begin{array}[]{c}Vr_{j}\\\
Wr_{j}\end{array}\right]=\left[\begin{array}[]{c}0\\\
y_{j}\end{array}\right],\quad\left[\begin{array}[]{cc}\hat{A}-\lambda_{j}I_{n}&\hat{B}\\\
\hat{C}&D\end{array}\right]\left[\begin{array}[]{c}r_{j}\\\
Wr_{j}\end{array}\right]=\left[\begin{array}[]{c}0\\\
y_{j}\end{array}\right],$

where $y_{j}:=(CV+DW)r_{j}$. This then implies that
$y_{j}=Z(\lambda_{j})Wr_{j}=\hat{Z}(\lambda_{j})Wr_{j}$, which together with
(3.10), (3.11) and $\hat{D}=D$ yields (3.9). ∎

Notice that Theorem 4 constructs a reduced-order system and a corresponding
certificate $\hat{X}$ for passivity, but not for strict passivity, since the
matrix ${\mathcal{W}}(\hat{X},\hat{\mathcal{M}})$ is positive semi-definite
and singular, while we would prefer to construct a lower order model with a
certificate for strict passivity.

##  4 Parameterized interpolants

In this section we combine the results of Sections 2 and 3 to propose a set of
parameterized interpolants that have the property that the interpolants have a
realization that is port-Hamiltonian and at the same time a passivity radius
that has a sufficiently large lower bound.

For this, we proceed as follows. Let $Z(s)$ be a strictly passive transfer
function of McMillan degree $n$, and suppose we are given a minimal model
${\mathcal{M}}:=\\{A,B,C,D\\}$ of $Z(s)$. We will then construct a lower order
model via the spectral zeros method explained in Section 3 but applied to a
so-called shifted transfer function :

$Z_{\xi}(s):=Z(s-\frac{\xi}{2})-\frac{\xi}{2}I_{m},\quad\mathrm{with\;model}\quad{\mathcal{M}}_{\xi}:=\\{A+\frac{\xi}{2}I_{n},B,C,D-\frac{\xi}{2}I_{m}\\},$
(4.1)

where $\xi$ is chosen in the open interval $(0,\Xi)$ of strictly passive
systems $Z_{\xi}(s)$ (see Theorem 3). We then solve the tangential
interpolation problem to produce a lower order model $\hat{Z}_{\xi}(s)$ of
degree $\hat{n}<n$ using
$\hat{Z}_{\xi}(\infty)=Z_{\xi}(\infty)=D-\frac{\xi}{2}I_{m}$ as well as
interpolation conditions on a subset of the spectral zeros of $Z_{\xi}(s)$ :

$Z_{\xi}(\sigma_{j})W_{\xi}r_{j}=\hat{Z}_{\xi}(\sigma_{j})W_{\xi}r_{j},\quad
r^{\mathsf{T}}_{j}W_{\xi}^{\mathsf{T}}Z_{\xi}(-\sigma_{j})=r^{\mathsf{T}}_{j}W_{\xi}^{\mathsf{T}}\hat{Z}_{\xi}(-\sigma_{j}),\quad
j=1,\ldots,\hat{n},$ (4.2)

where $(\sigma_{j},r_{j}),\;j=1,\ldots,\hat{n}$, are self-conjugate
(eigenvalue, eigenvector) pairs of the matrix $R_{\xi}$, chosen to have its
spectrum in the open right half plane, and which is obtained from the
deflating subspace equation

$\left[\begin{array}[]{cc|c}0&\\!A+\frac{\xi}{2}I_{n}-sI_{n}\\!&B\\\
\\!A^{\mathsf{T}}+\frac{\xi}{2}I_{n}+sI_{n}\\!&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D-\xi
I_{m}\end{array}\right]\left[\begin{array}[]{c}U_{\xi}\\\ V_{\xi}\\\ \hline\cr
W_{\xi}\end{array}\right]\\!=\\!\left[\begin{array}[]{c}V_{\xi}\\\
\\!-U_{\xi}\\\ \hline\cr 0\end{array}\right](R_{\xi}-sI_{\hat{n}}).$ (4.3)

It follows that these $2m\hat{n}+m^{2}$ real conditions completely define the
reduced-order model $\hat{Z}_{\xi}(s)$ and from Section 3 that a realization
of the reduced-order model is given by the quadruple

$\\{(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}(A+\frac{\xi}{2}I_{n})V_{\xi},(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}B,CV_{\xi},D-\frac{\xi}{2}I_{m}\\}$

which can also be written as

$\\{(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}AV_{\xi}+\frac{\xi}{2}I_{\hat{n}},(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}B,CV_{\xi},D-\frac{\xi}{2}I_{m}\\}.$

We now rewrite these conditions in terms of the original matrix $Z(s)$ and its
approximation $\hat{Z}(s)$ derived via this implicit shift technique.

######  Theorem 6.

Let ${\mathcal{M}}:=\\{A,B,C,D\\}$ be a minimal state-space realization of a
strictly passive transfer function $Z(s)$ of McMillan degree $n$, and let

$\Xi:=\sup_{\xi}\\{\xi\;|\;Z_{\xi}(s)\;\mathrm{is\;strictly\;passive}\\}.$

Then for any $\xi\in(0,\Xi)$, we consider an $\hat{n}$ dimensional deflating
subspace of the shifted pencil (4.3) corresponding to the spectrum of a real
matrix $R_{\xi}$ with eigenvalues in the right half plane. Then the matrices
$U_{\xi}$ and $V_{\xi}$ have full column rank $\hat{n}$, the matrix
$\hat{X}_{\xi}:=-U_{\xi}^{\mathsf{T}}V_{\xi}$ is symmetric and positive
definite, and the low order transfer function $\hat{Z}(s)$ with model
parameters

$\hat{\mathcal{M}}:=\\{\hat{A},\hat{B},\hat{C},\hat{D}\\}=\\{(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}AV_{\xi},(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}B,CV_{\xi},D\\}$
(4.4)

satisfies the interpolation conditions $Z(\infty)=\hat{Z}(\infty)=D$ and for
$\;j=1,...,\hat{n}$ :

$Z(\sigma_{j}-\xi/2)W_{\xi}r_{j}=\hat{Z}(\sigma_{j}-\xi/2)W_{\xi}r_{j},\quad
r^{\mathsf{T}}_{j}W_{\xi}^{\mathsf{T}}Z(-\sigma_{j}-\xi/2)=r^{\mathsf{T}}_{j}W_{\xi}^{\mathsf{T}}\hat{Z}(-\sigma_{j}-\xi/2),$
(4.5)

where $(\sigma_{j},r_{j}),\;j=1,...,\hat{n}$, are self-conjugate (eigenvalue,
eigenvector) pairs of the matrix $R_{\xi}$. Moreover, the matrix
$\hat{X}_{\xi}$ is a certificate for the LMI

${\mathcal{W}}(\hat{X}_{\xi},\hat{\mathcal{M}}):=\left[\begin{array}[]{cc}-\hat{X}_{\xi}\hat{A}-\hat{A}^{\mathsf{T}}\hat{X}_{\xi}&\hat{C}^{\mathsf{T}}-\hat{X}_{\xi}\hat{B}\\\
\hat{C}-\hat{B}^{\mathsf{T}}\hat{X}_{\xi}&\hat{D}^{\mathsf{T}}+\hat{D}\end{array}\right]\succeq\xi.\operatorname{diag}(\hat{X}_{\xi},I_{m})\succ
0$

and $\xi/2$ is a lower bound for the passivity radius of the normalized port-
Hamiltonian realization
$\hat{\mathcal{M}}_{T_{\xi}}:=\\{\hat{J}-\hat{R},\hat{G}-\hat{P},(\hat{G}+\hat{P})^{\mathsf{T}},\hat{N}+\hat{S}\\}$
obtained using $\hat{X}_{\xi}=T_{\xi}^{\mathsf{T}}T_{\xi}$ via the state-space
transformation

$\left[\begin{array}[]{cccc}T_{\xi}&0\\\
0&I_{m}\end{array}\right]\left[\begin{array}[]{cccc}\hat{A}&\hat{B}\\\
\hat{C}&\hat{D}\end{array}\right]\left[\begin{array}[]{cccc}T_{\xi}^{-1}&0\\\
0&I_{m}\end{array}\right]=\left[\begin{array}[]{cc}\hat{J}-\hat{R}&\hat{G}-\hat{P}\\\
(\hat{G}+\hat{P})^{\mathsf{T}}&\hat{N}+\hat{S}\end{array}\right].$ (4.6)

###### Proof.

It follows from the strict passivity of $Z_{\xi}(s)$ for any $\xi$ in the open
interval $(0,\Xi)$ that $\hat{Z}_{\xi}(s)$ constructed using (4.3) and (4.4),
satisfies the conditions of Theorem 4 and Remark 2. Therefore, the matrix
$\hat{X}_{\xi}:=-U_{\xi}^{\mathsf{T}}V_{\xi}$ is symmetric and positive
definite. It then follows that the projected system

$\hat{\mathcal{M}}_{\xi}:=\\{(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}A_{\xi}V_{\xi},(U_{\xi}^{\mathsf{T}}V_{\xi})^{-1}U_{\xi}^{\mathsf{T}}B_{\xi},C_{\xi}V_{\xi},D_{\xi}\\}$

satisfies Theorem 4 with $\hat{X}_{\xi}:=-U_{\xi}^{\mathsf{T}}V_{\xi}$ and
hence we have

${\mathcal{W}}(\hat{X}_{\xi},\hat{\mathcal{M}}_{\xi})=\left[\begin{array}[]{cc}-\hat{X}_{\xi}\hat{A}_{\xi}-\hat{A}^{\mathsf{T}}_{\xi}\hat{X}_{\xi}&\hat{C}^{\mathsf{T}}_{\xi}-\hat{X}_{\xi}\hat{B}_{\xi}\\\
\hat{C}_{\xi}-\hat{B}^{\mathsf{T}}_{\xi}\hat{X}_{\xi}&\hat{D}_{\xi}+\hat{D}^{\mathsf{T}}_{\xi}\end{array}\right]\succeq
0.$

By using the relations between
$\hat{\mathcal{M}}_{\xi}=\\{\hat{A}_{\xi},\hat{B}_{\xi},\hat{C}_{\xi},\hat{D}_{\xi}\\}=\\{\hat{A}+\frac{\xi}{2}I_{n},\hat{B},\hat{C},\hat{D}-\frac{\xi}{2}I_{m}\\}$
and $\hat{\mathcal{M}}=\\{\hat{A},\hat{B},\hat{C},\hat{D}\\}$, we obtain the
LMI

${\mathcal{W}}(\hat{X}_{\xi},\hat{\mathcal{M}})=\left[\begin{array}[]{cc}-\hat{X}_{\xi}\hat{A}-\hat{A}^{\mathsf{T}}\hat{X}_{\xi}&\hat{C}^{\mathsf{T}}-\hat{X}_{\xi}\hat{B}\\\
\hat{C}-\hat{B}^{\mathsf{T}}\hat{X}_{\xi}&\hat{D}+\hat{D}^{\mathsf{T}}\end{array}\right]\succeq\xi\left[\begin{array}[]{cc}\hat{X}_{\xi}&0\\\
0&I_{m}\end{array}\right]\succ 0$ (4.7)

which implies that the transformed port-Hamiltonian system (4.6) has a
passivity radius at least as large as $\frac{\xi}{2}$ since it follows from
(4.7) and $\hat{X}_{\xi}=T_{\xi}^{\mathsf{T}}T_{\xi}$, that

$\frac{1}{2}\left[\begin{array}[]{cccc}T_{\xi}^{\mathsf{-T}}&0\\\
0&I_{m}\end{array}\right]{\mathcal{W}}(\hat{X}_{\xi},\hat{\mathcal{M}})\left[\begin{array}[]{cccc}T_{\xi}^{-1}&0\\\
0&I_{m}\end{array}\right]=\left[\begin{array}[]{cc}\hat{R}&\hat{P}\\\
\hat{P}^{\mathsf{T}}&\hat{S}\end{array}\right]\succ\frac{\xi}{2}I_{\hat{n}+m}.$

The translation of interpolation conditions on the shifted system towards
similar conditions on the original system follows directly from the identity
(4.1). ∎

It follows from the above theorem and from Theorem 2 that in order to have an
optimal passivity radius for the reduced-order model, one should choose to put
it in the normalized port-Hamiltonian form
$\\{T_{\xi}\hat{A}T_{\xi}^{-1},T_{\xi}\hat{B},\hat{C}T_{\xi}^{-1},\hat{D}\\}$.

######  Remark 3.

It follows from Theorem 4 that when choosing
$\Xi<\xi<\lambda_{\min}(D^{\mathsf{T}}+D)$, the pencil $S_{\xi}(s)$ may still
have a deflating subspace (4.3) where $\hat{X}_{\xi}$ is positive definite,
and hence yield a strictly passive reduced-order model. If this is the case,
we will be able to increase the passivity radius even further. This
flexibility will be used in the section on numerical examples.

######  Remark 4.

Notice that the interpolation points $\\{\sigma_{j},j=1,\ldots,\hat{n}\\}$ and
$\\{-\sigma_{j},j=1,\ldots,\hat{n}\\}$ of the shifted system
$\hat{Z}_{\xi}(s)$ are mirror images of each other with respect to the origin,
but this is not true anymore for the interpolation points
$\\{\sigma_{j}-\xi/2,j=1,\ldots,\hat{n}\\}$ and
$\\{-\sigma_{j}-\xi/2,j=1,\ldots,\hat{n}\\}$ of the original system
$\hat{Z}(s)$. Moreover, since the interpolation points
$\\{\sigma_{j},j=1,\ldots,\hat{n}\\}$ are still in the open right half plane,
the shifted interpolation conditions have the tendency to approximate better
the transfer function in the left half plane.

##  5 Choosing the spectral zeros

In this section we look at the selection of zeros and the effect of (near)
non-minimality of the transfer function. If we want to select particular
spectral zeros, it is convenient to compute the individual corresponding
eigenvectors :

$\left[\begin{array}[]{cc|c}0&A-\lambda_{j}I_{n}&B\\\
A^{\mathsf{T}}+\lambda_{j}I_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right]\left[\begin{array}[]{c}Ur_{j}\\\
Vr_{j}\\\ \hline\cr Wr_{j}\end{array}\right]=0,\quad Rr_{j}=\lambda_{j}r_{j}.$
(5.1)

It follows from the proof of Theorem 5 that the interpolation condition
becomes

$y_{j}:=\left[CVr_{j}+DWr_{j}\right]=Z(\lambda_{j})Wr_{j}=\hat{Z}(\lambda_{j})Wr_{j}.$
(5.2)

Since $y_{j}$ and $Wr_{j}$ are both bounded quantities, $\lambda_{j}$ can not
be a pole of $Z(s)$ unless it is also a decoupling zero, implying that the
system is not minimal. More formally, let $\cal N$ be the unobservable
subspace of the pair $(A,C)$ then

$A{\cal N}\subset{\cal N},\quad C{\cal N}=0.$

This implies that $0\oplus{\cal N}\oplus 0$ is a deflating subspace of $S(s)$
with as spectrum the unobservable modes of the pair $(A,C)$. Choosing a vector
in that deflating subspace yields $Ur_{j}=0$ and hence also a singular matrix
$\hat{X}$. Moreover, one then has $Wr_{j}=0$ and the interpolation condition
(5.2) then vanishes. A similar reasoning on the dual system implies that the
same problem occurs when using an uncontrollable mode of the pair $(A,B)$.
Therefore it is recommended to stay away from nearly uncontrollable or
unobservable modes when selecting spectral zeros as interpolation points. If
we make sure that $\hat{X}$ has large eigenvalues, then we will stay away from
non-minimality in the reduced-order model, and the interpolation conditions
$\eqref{interpolate}$ will be well defined. Moreover, it makes sense to
“maximize” $\hat{X}$ since it is the Hamiltonian storage function of the
projected system: maximizing $\hat{X}$ can indeed be viewed as finding the
dominant restriction of the Hamiltonian $X$.

We used the following procedure to construct a nearly optimal selection of
interpolating spectral zeros. Assume that we computed the full matrix
$X:=-U^{\mathsf{T}}V$, which is symmetric and positive definite. If we perform
the Cholesky decomposition with pivoting on this matrix, then the leading
$\hat{n}\times\hat{n}$ submatrix $\hat{X}$ corresponds to a subset of
$\hat{n}$ eigenvectors that is nearly optimal (a truly optimal selection would
require to verify all possible symmetric permutations). In practice this
“greedy” ordering of the spectral zeros works reasonably well on the examples
we tried. We should point out that if one desires a real lower order model,
then the pivoting strategy should also make sure that the selected spectral
zeros form a self-conjugate set, but that is easy to obtain via a post-
processing of the greedy ordering : it amounts to looking for a leading subset
of $\hat{n}$ self-conjugate spectral zeros in the preliminary ordered complex
zeros. We also point out that this selection procedure can also be implemented
on a partial set of computed eigenvectors and spectral zeros, such as those
one would compute using a Krylov-Schur method for large-scale problems (see
e.g. [7]) combined with implicit filtering of undesired spectral zeros. Such
large-scale issues, though, are beyond the scope of this paper.

##  6 Using the robustness property

It follows from Section 4 that it is indicated to choose $\xi\in(0,\Xi)$ as
large as possible, since this will yield an interpolant with a certificate for
a larger passivity radius. This means that in that coordinate system we can
allow for larger perturbations and still preserve passivity of the reduced-
order model. We can therefore expect to have more freedom in the numerical
implementation of any algorithm computing the deflating subspace described in
Theorem 6 or on the flexibility of its stopping criterion.

We first show that for a strictly passive system, there are many possibilities
to construct strictly passive lower order models and that the corresponding
projectors form an open set.

######  Theorem 7.

Let ${\mathcal{M}}:=\\{A,B,C,D\\}$ be a minimal state-space model for a
strictly passive transfer function $Z(s)$ of McMillan degree $n$. Let $X\succ
0$ be a certificate for the LMI that ensures that $Z(s)$ is strictly passive :

${\mathcal{W}}(X,{\mathcal{M}})=\left[\begin{array}[]{cc}-X&0\\\
0&I_{m}\end{array}\right]\left[\begin{array}[]{cc}A&B\\\
C&D\end{array}\right]+\left[\begin{array}[]{cc}A^{\mathsf{T}}&C^{\mathsf{T}}\\\
B^{\mathsf{T}}&D^{\mathsf{T}}\end{array}\right]\left[\begin{array}[]{cc}-X&0\\\
0&I_{m}\end{array}\right]\succ 0.$ (6.1)

If we choose any matrix $V\in\mathbb{R}^{n\times\hat{n}}$ of full column rank
$\hat{n}$, and compute $U:=-XV\hat{X}^{-1}$, where
$\hat{X}:=V^{\mathsf{T}}XV$, then $U^{\mathsf{T}}V=-I_{\hat{n}}$ and the
system

$\hat{\mathcal{M}}:=\\{\hat{A},\hat{B},\hat{C},\hat{D}\\}=\\{(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}AV,(U^{\mathsf{T}}V)^{-1}U^{\mathsf{T}}B,CV,D\\}$
(6.2)

is a strictly passive lower order model of degree $\hat{n}$.

###### Proof.

It follows from (6.1) that

$\left[\begin{array}[]{cc}\\!-V^{\mathsf{T}}X\\!&0\\\
0&\\!I_{m}\\!\end{array}\right]\left[\begin{array}[]{cc}A&B\\\
C&D\end{array}\right]\left[\begin{array}[]{cc}V&0\\\
0&\\!I_{m}\\!\end{array}\right]+\left[\begin{array}[]{cc}\\!V^{\mathsf{T}}\\!&0\\\
0&I_{m}\end{array}\right]\left[\begin{array}[]{cc}A^{\mathsf{T}}&C^{\mathsf{T}}\\\
B^{\mathsf{T}}&D^{\mathsf{T}}\end{array}\right]\left[\begin{array}[]{cc}\\!-XV\\!&0\\\
0&\\!I_{m}\\!\end{array}\right]\succ 0.$ (6.3)

Using $\hat{X}=V^{\mathsf{T}}XV$, $U\hat{X}=-XV$ and
$U^{\mathsf{T}}V=-I_{\hat{n}}$, we can rewrite this as

$\left[\begin{array}[]{cc}-\hat{X}&0\\\
0&I_{m}\end{array}\right]\left[\begin{array}[]{cc}\hat{A}&\hat{B}\\\
\hat{C}&\hat{D}\end{array}\right]+\left[\begin{array}[]{cc}\hat{A}^{\mathsf{T}}&\hat{C}^{\mathsf{T}}\\\
\hat{B}^{\mathsf{T}}&\hat{D}^{\mathsf{T}}\end{array}\right]\left[\begin{array}[]{cc}-\hat{X}&0\\\
0&I_{m}\end{array}\right]\succ 0$

which proves the strict passivity of the lower order model, since
$\hat{X}\succ 0$. Moreover, the matrices $U$, $V$, $\hat{X}$ and
$U^{\mathsf{T}}V$ have full rank $\hat{n}$ by construction and this is
maintained in an open neighborhood of $U$ and $V$. Therefore the matrix
inequality (6.3) is still valid and the constructed reduced-order models in a
sufficiently small neighborhood of (6.2) are strictly passive. ∎

Let us suppose now that the deflating subspace described in (4.3) was
inaccurate, either due to roundoff, or due to early termination of an
iterative process to compute it. If we denote the computed quantities as
$\widetilde{U}$, $\widetilde{V}$ and $\widetilde{W}$, then we can construct
$\widetilde{R}$ and residuals $\Delta_{U}$, $\Delta_{V}$ and $\Delta_{W}$ such
that the following equation holds

$\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\\
A^{\mathsf{T}}+sI_{n}&0&C^{\mathsf{T}}\\\ \hline\cr
B^{\mathsf{T}}&C&D^{\mathsf{T}}+D\end{array}\right]\left[\begin{array}[]{c}\widetilde{U}\\\
\widetilde{V}\\\
\hline\cr\widetilde{W}\end{array}\right]=\left[\begin{array}[]{c}\widetilde{V}\\\
-\widetilde{U}\\\ \hline\cr
0\end{array}\right](\widetilde{R}-sI_{\hat{n}})+\left[\begin{array}[]{c}\Delta_{U}\\\
\Delta_{V}\\\ \hline\cr\Delta_{W}\end{array}\right].$ (6.4)

Let us also denote the computed projected system as

$\widetilde{\mathcal{M}}:=\\{\widetilde{A},\widetilde{B},\widetilde{C},\widetilde{D}\\}:=\\{(\widetilde{U}^{\mathsf{T}}\widetilde{V})^{-1}\widetilde{U}^{\mathsf{T}}A\widetilde{V},(\widetilde{U}^{\mathsf{T}}\widetilde{V})^{-1}\widetilde{U}^{\mathsf{T}}B,C\widetilde{V},D\\}.$

If we define $\widetilde{X}:=-\widetilde{U}^{\mathsf{T}}\widetilde{V}$, then
if follows from these equations that

$(\widetilde{X}^{\mathsf{T}}-\widetilde{X})\widetilde{R}+\widetilde{R}^{\mathsf{T}}(\widetilde{X}^{\mathsf{T}}-\widetilde{X})=(\Delta_{U}^{\mathsf{T}}\widetilde{U}+\Delta_{V}^{\mathsf{T}}\widetilde{V}+\Delta_{W}^{\mathsf{T}}\widetilde{W})-(\widetilde{U}^{\mathsf{T}}\Delta_{U}+\widetilde{V}^{\mathsf{T}}\Delta_{V}+\widetilde{W}^{\mathsf{T}}\Delta_{W}),$
(6.5)

which implies that $\widetilde{X}$ is nearly symmetric, and that the following
matrix is nearly positive definite :

$\left[\begin{array}[]{cc}-\widetilde{X}\widetilde{A}-\widetilde{A}^{\mathsf{T}}\widetilde{X}^{\mathsf{T}}&\widetilde{C}^{\mathsf{T}}-\widetilde{X}\widetilde{B}\\\
\widetilde{C}-\widetilde{B}^{\mathsf{T}}\widetilde{X}^{\mathsf{T}}&\widetilde{D}^{\mathsf{T}}+\widetilde{D}\end{array}\right]=\left[\begin{array}[]{cc}\widetilde{U}^{\mathsf{T}}A\widetilde{V}+\widetilde{V}^{\mathsf{T}}A^{\mathsf{T}}\widetilde{U}&\widetilde{V}^{\mathsf{T}}C^{\mathsf{T}}+\widetilde{U}^{\mathsf{T}}B\\\
C\widetilde{V}+B^{\mathsf{T}}\widetilde{U}&D^{\mathsf{T}}+D\end{array}\right]$
$=\left[\begin{array}[]{cc}\widetilde{W}^{\mathsf{T}}\\\
-I_{m}\end{array}\right](\widetilde{D}^{\mathsf{T}}+\widetilde{D})\left[\begin{array}[]{cc}\widetilde{W}&-I_{m}\end{array}\right]+\left[\begin{array}[]{cc}\Delta&\Delta_{W}\\\
\Delta_{W}^{\mathsf{T}}&0\end{array}\right]\succeq 0,$

where
$\Delta=(\widetilde{X}^{\mathsf{T}}-\widetilde{X})\widetilde{R}+\widetilde{U}^{\mathsf{T}}\Delta_{U}+\widetilde{V}^{\mathsf{T}}\Delta_{V}-\Delta_{W}^{\mathsf{T}}\widetilde{W}$
is symmetric, because of (6.5).

Notice that this is not a valid passivity LMI since $\widetilde{X}$ is not
symmetric. But if we replace $\widetilde{X}$ by its symmetric part
$\widetilde{X}_{s}=\frac{1}{2}(\widetilde{X}+\widetilde{X}^{\mathsf{T}})$ then
we obtain, using
$\widetilde{X}_{a}=\frac{1}{2}(\widetilde{X}-\widetilde{X}^{\mathsf{T}})$

$\left[\begin{array}[]{cc}\\!-\widetilde{X}_{s}\widetilde{A}-\widetilde{A}^{\mathsf{T}}\widetilde{X}_{s}\\!&\\!\widetilde{C}^{\mathsf{T}}-\widetilde{X}_{s}\widetilde{B}\\!\\\
\\!\widetilde{C}-\widetilde{B}^{\mathsf{T}}\widetilde{X}_{s}\\!&\widetilde{D}^{\mathsf{T}}+\widetilde{D}\end{array}\right]\\!=\\!\left[\begin{array}[]{cc}\widetilde{W}^{\mathsf{T}}\\\
\\!-I_{m}\\!\end{array}\right](\widetilde{D}^{\mathsf{T}}+\widetilde{D})\left[\begin{array}[]{cc}\\!\widetilde{W}\\!&\\!-I_{m}\\!\end{array}\right]+\left[\begin{array}[]{cc}\Delta_{11}&\\!\Delta_{12}\\!\\\
\Delta_{12}^{\mathsf{T}}&0\end{array}\right]\succeq 0,$

where
$\Delta_{11}=\Delta-\widetilde{X}_{a}\widetilde{A}-\widetilde{A}^{\mathsf{T}}\widetilde{X}_{a}$
and $\Delta_{12}=\Delta_{W}-\widetilde{X}_{a}\widetilde{B}$. Notice that
$\widetilde{X}_{a}$ is a solution of the Lyapunov-like equation (6.5) and
hence that the perturbation of the above passivity LMI is of the order of the
residual in (6.4).

This shows that if we have a robustness margin in the unperturbed system, in
the sense that its passivity radius is bounded away from 0, then strict
passivity is maintained for a reasonably large residual in (6). We can thus
apply these ideas to the technique of shifted interpolation and guarantee that
the perturbations induced by the numerical algorithm do not destroy the strict
passivity of the projected model. Notice that when using iterative algorithms
for large-scale problems, such robustness properties may come in handy since
we may allow for early termination of iterative schemes, provided the
resulting perturbation lies within the robustness bounds.

##  7 Numerical experiments

In this section, we illustrate the proposed methodology to construct passive
reduced-order models by means of two numerical examples. All the experiments
were conducted using MATLAB®2020b.

###  7.1 RLC circuit:

We first illustrate the results of the parameterized interpolation technique
by applying it to the $200^{th}$ order single-input/single-output model of a
circuit described in [12], where $100$ electrical capacitances, inductors, and
resistances are interconnected. The limiting value $\Xi\approx 0.56$ for the
parameter $\xi$ was estimated using a mesh of equidistant points in the
interval $[0,\Xi_{ub}]$, where $\Xi_{ub}$ is a conservative upper bound
computed from the spectrum of $A$ (see [13]). We applied the selection
procedure of spectral zeros described in Section 5 for lower order degrees
$k=\\{2,4,\ldots,20\\}$ and for equidistant shifts $\xi\in[0,\Xi]$.

Figure 7.1: Degree $6$ approximation of 200th order RLC network. Top left:
spectral zeros (dots) and right half plane interpolation sets of points (+).
Black and green indicate the interpolation points corresponding to $\xi=0$ and
$\xi=\Xi$, respectively. Top right: relative $H_{\infty}$-error norm as a
function of $\xi$. Bottom left: singular value plot of the original system
(SYS) and its best approximation (sys). Bottom right: singular value plot of
the corresponding error system.

In Figure 7.1, we give the results of the low order model of degree $6$ for
equidistant shifts $\xi\in[0,0.56]$. The top-left plot shows the original
spectral zeros (in blue dots) and the selected right half plane interpolation
points in magenta color for different values of $\xi$. Moreover, low-intensity
magenta color $`+`$ belongs to lower values of $\xi$; likewise, high-intensity
color belongs to larger values of $\xi$. One can see that the interpolation
points are close to the original spectral zeros but with a shift towards the
imaginary axis as $\xi$ increases. The top right plot gives the relative
$H_{\infty}$-error norm as a function of $\xi$. One can see that the errors
depend in a non-smooth manner on the parameter $\xi$, which is not so
surprising since the selected interpolation points also depend on $\xi$. It is
to be noted, though, that there is a general decreasing trend of the relative
error as a function of $\xi$. This is also the case for the other low-order
models we constructed. The bottom two plots give the singular value plot of
the original system (SYS) and its best approximation (sys), and the singular
value plot of the corresponding error system (SYS-sys), respectively.

Figure 7.2: Bounds for the low order error norms. The solid lines lines give
upper and lower bounds for the models constructed for different values of
$\xi$. The dashed lines give provable upper and lower bounds for the optimal
$H_{\infty}$ reduced-order model.

In order to show the effect of the order selection, we give in Figure 7.2 a
plot of lower and upper values of the achieved relative $H_{\infty}$-errors
for the different values of $\xi$, as a function of the order $k$. In
comparison, we also included provable upper and lower bounds of the relative
error for the optimal $H_{\infty}$ approximation of the given system. Note
that the upper bound is computed using the standard balanced truncation
method, see, e.g., [2], and the lower error bound for reduced models of order
$r$ can be determined by $\sigma_{k+1}$, where $\sigma_{k+1}$ is the $(k+1)$th
largest singular value of the original system [17]. It is clear from this plot
that our selection procedure is far from optimal. One should be aware, though,
that our procedure is restricted to lower order systems that are passive and
are generated by interpolation of special sets of points, which is a
restrictive constraint.

###  7.2 Random example:

The second example is a random port-Hamiltonian system
${\mathcal{M}}:=\\{J-R,G-P,G^{\mathsf{T}}+P^{\mathsf{T}},N+S\\}$ with
$X$-passivity radius $\rho_{\mathcal{M}}(I_{n})=0.5$. The state-space model
has state dimension $n=6$ and input/output dimension $m=2$. It was generated
by constructing a random symmetric matrix
${\mathcal{W}}:=\left[\begin{smallmatrix}R&P\\\
P^{\mathsf{T}}&S\end{smallmatrix}\right]$ with smallest eigenvalue
$\lambda_{\min}=0.5$ and a random anti-symmetric matrix
${\mathcal{V}}:=\left[\begin{smallmatrix}-J&-G\\\
G^{\mathsf{T}}&N\end{smallmatrix}\right]$. We then applied equidistant shifts
$\xi\in[0,\Xi]$ and computed reduced-order models of degree $\hat{n}=4$, based
on the parameterized method explained in Theorem 6. But based on Remark 3, we
also took values of the shift $\xi>\Xi$, for as long as the construction of a
positive definite matrix $\hat{X}_{\xi}$ was possible (which implies that
$\hat{A}_{\xi}$ is still stable and
$\hat{D}^{\mathsf{T}}_{\xi}+\hat{D}_{\xi}\succ 0$).

In Figure 7.3, the top-left plot shows the original spectral zeros (in blue
dots) and selected right half plane interpolation points as magenta-colored
$`+`$. One can see that as a function of $\xi$, the choice of four
interpolation points is now much closer to each other than in the previous
example. The top right plot gives the relative $H_{\infty}$-error norm as a
function of $\xi$, and there also, one observes a smoother behavior since
essentially the same interpolation points are being used. We point out here
that in this plot, the blue crosses correspond to the values of
$\xi\in[0,\Xi]$, whereas the red circles, correspond to the values of
$\xi>\Xi$. The bottom two plots give the singular value plots of the original
system (SYS) and its best approximation (sys) and the singular value plots of
the corresponding error system (SYS-sys), respectively.

Figure 7.3: Degree $4$ approximation of random $6$th order network. Top left:
spectral zeros (dots) and interpolation sets of points (crosses). Black and
green indicate the interpolation points corresponding to $\xi=0$ and
$\xi\approx 1.60\Xi$, respectively. Top right: relative $H_{\infty}$-error
norm as a function of $\xi$. Bottom left: singular value plots of the original
system (SYS) and its best approximation (sys). Bottom right: singular value
plots of the corresponding error system (SYS-sys).

##  8 Concluding remarks

In this paper we developed a parameterized model reduction method based on the
interpolation of the transfer function $Z(s)$ in a subset of the so-called
spectral zeros of $Z(s)$. The parameterization lies in the fact that we
consider now spectral zeros of shifted systems $Z_{\xi}(s)$, rather than the
original transfer function $Z(s)$. Although the method is theoretically based
on interpolation techniques, the algorithm itself is based on the computation
of particular deflating subspaces of the “Hamiltonian” pencils associated with
the shifted models $Z_{\xi}(s)$. It was also shown that these deflating
subspaces do not need to be computed exactly, since the bounds on the
passivity radius of the projected systems, gives a certain flexibility in the
accuracy needed for the eigenspace computation. In this paper, we also
proposed a new procedure for the selection of spectral zeros used as
interpolation points for a lower order model that is a good approximation.

## Acknowledgment

The research was performed during two visits of the third author to the Max
Planck Institute in Magdeburg.

## References

  * [1] B.D.O. Anderson and S. Vongpanitlerd,  Network Analysis and Synthesis.  Dover Pblications, 2006.
  * [2] A.C. Antoulas,  Approximation of Large-Scale Dynamical Systems.  Advances in Design and Control, Vol. DC-06, SIAM, Philadelphia, 2005.
  * [3] A.C. Antoulas,  A new result on passivity preserving model reduction.  Systems and Control Letters, Vol. 54:361-374, 2005.
  * [4] C. Beattie, V. Mehrmann and H. Xu,  Port-Hamiltonian realizations of linear time-invariant systems.  arXiv:2201.05355, 2022.
  * [5] P. Benner, V. Mehrmann, and D.C. Sorensen (Eds.),  Dimension Reduction of Large-Scale Systems.  Lecture Notes in Computational Science and Engineering, Vol. 45, Springer-Verlag, Berlin/Heidelberg, June 2005.
  * [6] P. Benner, A. Cohen, M. Ohlberger, and K. Willcox (Eds.),  Model Reduction and Approximation: Theory and Algorithms.  Computational Science and Engineering, Vol. 15, SIAM Publications, Philadelphia, PA, 2017.
  * [7] P. Benner, H. Fassbender and M. Stoll,  A Hamiltonian Krylov–Schur-type method based on the symplectic Lanczos process.  Linear Algebra and its Applications, Vol. 435:578-600, 2011.
  * [8] P. Benner, and H. Fassbender,  Numerical methods for passivity preserving model reduction.  at-Automatisierungstechnik, Vol. 54(4):153-160, 2006 (in German).
  * [9] G. Fanizza, J. Karlsson, A. Lindquist, and R. Nagamune,  Passivity-preserving model reduction by analytic interpolation.  Linear Alg. Appl., Vol. 425:608-633, 2007.
  * [10] K. Gallivan, A. Vandendorpe, P. Van Dooren,  Model reduction of MIMO systems via tangential interpolation.  SIAM J. Matrix Anal. Appl., Vol. 26:328-349, 2004.
  * [11] Y. Genin, Y. Hanchez, Y. Nesterov, R. Stefan, P. Van Dooren and S. Xu,  Positivity and linear matrix inequalities.  European Journal of Control, Vol. 8(3):275-298, 2002.
  * [12] S. Gugercin, A.C. Antoulas,  A survey of balancing methods for model reduction.  in: Proc. European Control Conf. ECC 2003, Cambridge, UK, 2003, CD Rom.
  * [13] V. Mehrmann and P. Van Dooren,  Optimal robustness of port-Hamiltonian systems.  SIAM J. Matr. Anal. Appl., Vol. 41(1):134-151, 2019.
  * [14] A.C. Sorensen,  Passivity preserving model reduction via interpolation of spectral zeros.  Systems and Control Letters, Vol. 54:347-360, 2005.
  * [15] A.J. van der Schaft and D. Jeltsema,  Port-Hamiltonian Systems Theory, an Introductory Survey.  Foundations and Trends in Systems and Control, Vol. 1:173-378, 2014.
  * [16] J. C. Willems,  Dissipative dynamical systems – Part II: Linear systems with quadratic supply rates.  Arch. Ration. Mech. Anal., Vol. 45:352-393, 1972.
  * [17] K. Glover,  All optimal Hankel-norm approximations of linear multivariable systems and their L∞-error norms.  Internat. J. Control, Vol. 39:1115-1193, 1984.