N Le Neindre

R Bougault

S Barlini

E Bonnet

B Borderie

G Casini

A Chbihi

P Edelbruck

J D Frankland

D Gruyer

E Legouée

O Lopez

P Marini

M Pârlog

G Pasquali

M Petcu

M F Rivet

F Salomon

E Vient

R Alba

G Baiocco

L Bardelli

M Bini

R Borcea

M Bruno

S Carboni

M Cinausero

I Cruceru

M Degerlier

J A Dueñas

K Gasior

F Gramegna

A Grzeszczuk

M Kamuda

T Kozik

V Kravchuk

I Lombardo

C Maiolino

T Marchi

L Morelli

F Negoita

A Olmi

H Petrascu

S Piantelli

G Poggi

E Rosato

D Santonocito

G Spadaccini

A A Stefanini

T Twaróg

M Vigilante

Comparison of charged particle identification using pulse shape discrimination and ∆E-E methods between front and rear side injection in silicon detectors

Keywords: Pulse shape discrimination, ∆E-E, Silicon detector, Charge signal, Current signal PACS: 29.30.Ep, 29.30.-h, 29.40.-n, 29.40.Wk, 84.30.Sk

Comparison of charged particle identification using pulse shape discrimination and ∆E-E methods between front and rear side injection in silicon detectors

Introduction

The goal for the next generation of nuclear physics experiments will be to take maximum advantage of the increasing number of radioactive heavy-ion beams delivered by the existing and planned accelerator facilities. The nuclear structure as well as the thermodynamics and dynamics of excited exotic nuclei far from the stability valley can be studied. In particular, for the equation of state of nuclear matter it will be possible to explore the temperature and density dependence of the symmetry energy term on isospin. The aim of the FAZIA collaboration [START_REF]Web site Fazia collaboration[END_REF] is the definition and construction of a detector system for this purpose. A high granularity and large acceptance 4π apparatus for a detailed study of heavy-ion collisions in a wide range of bombarding energies (from a few to about 50 A MeV) is planned. In order to study the isospin degree of freedom, it is very important to extend as far as possible the isotopic identification capabilities of the detectors towards heavy fragments while keeping the energy thresholds very low. With this objective, few years ago the FAZIA collaboration started an R&D program, based on the optimization of the performances of three-layers Si-Si-CsI telescopes and sizeable improvements have been obtained for both the standard ∆E-E technique and the Pulse Shape Analysis [START_REF] Pausch | [END_REF] (hereafter PSA) in each of the silicon stages. Most efforts have been spent on the latter for ions which are stopped in the first Si-layer and cannot be identified via the ∆E-E method. Indeed, because of the different stopping powers, different particles with the same kinetic energy, or the same velocity, produce different energy loss profiles along the detector depth and this results in different charge collection times, i.e. in different pulse shapes [3]. Both techniques, and in particular the PSA, have benefited from the digital signal processing made possible by a new digital sampling electronics [4], [5], [6]. Other groups [START_REF] Pausch | [END_REF], [7], [8], [9], [10], [START_REF] Amorini | IEEE Nuclear Science Symposium Conference Record N[END_REF], [START_REF] Mutterer | [END_REF], [13], [14], [15], [16] have also investigated the behaviour of silicon detectors for PSA applications in previous years.

The aim of the present paper is to compare PSA techniques for two different particle injection. Indeed, the shape of the signals produced by particles crossing a semiconductor detector is very sensitive to the strength and the configuration of the electric field F(x) inside the detector. Therefore, as expected, we have observed a different behaviour when nuclear products enter the silicon np junction through the front side with the high electric field or through the rear side with the low electric field [17]. For a specified ion, the shape (as a function of time) of the induced signal depends, besides the specific rate of energy loss and the electric field, on the detector capacitance and preamplifier characteristics. Using exactly the same detector and preamplifier for particle incidence on the front or on the rear side should provide a better understanding of the signal shapes [17]. In literature both configurations have been successfully tested [START_REF] Pausch | [END_REF], [6], [9], but even when both mountings were used (e.g. in [START_REF] Pausch | [END_REF]) they have not been compared for the same detector and electronics. The front configuration presents a priori a big advantage if one wants to perform, besides the pulse shape, a Time of Flight (ToF) measurement, since the signals are faster. But the gain in ToF identification can imply a spoiling of PSA and a compromise should be found. In this experiment for the first time the very same detectors, preamplifiers, applied voltages and electronics have been used for both settings, allowing a fair comparison exactly in the same conditions, including reaction product distributions. PSA correlations and their qualities in terms of energy thresholds and shapes will be investigated.

Experimental setup

A 35 A MeV 84 Kr beam, delivered by the "Ciclotrone Superconduttore"(CS) of the "Laboratori Nazionali del Sud" (LNS) in Catania, was used for bombarding thin targets of 124 Sn, 120 Sn, 112 Sn, and nat Au. All results presented in this paper refer to the reaction products obtained by adding up the data of those reactions, in order to obtain larger statistics.

Telescope mounting

In this experiment we used the same telescope configuration as the one foreseen for the future FAZIA 4π array. The telescope composed of Si 300 µm-Si 500 µm-CsI(Tl) 10 cm, was mounted in the "Ciclope" scattering chamber at LNS, at a distance from target of ≃1 m and at ≃5.4 o polar angle with respect to the beam line, slightly beyond the grazing angles, i.e. the scattering angles corresponding to very peripheral projectile-target nuclear collisions, where the reactions mechanism concentrates the largest variety of products. Conversely to a previous experiment [6] a CsI detector was added to reject particles punching through the second silicon Si2. Some characteristics of the telescope and of its dedicated electronics are given in Table 1.

Properties of the detectors and their electronics

The silicon detectors come from random cut wafers. They are ion-implanted, of the neutron transmutation doped (n-TD) type, with bulk resistivity values in the range 3000-4000 Ωcm and an active area of 20×20 mm 2 . They were manufactured by FBK (Trento, Italy) in transmission mounting, with dead layers on both sides of ≃500-800 nm to improve their robustness. The 10 cm long CsI crystal is produced by Amcrys. Its lateral surfaces have been wrapped with the new Vikuiti T M ESR 3M reflecting polymer, while its back square face has been glued to a custom design photodiode produced by FBK.

Each element of the FAZIA telescopes (silicon detectors and CsI scintillators) was coupled to custom-built Front End Electronics (FEE), consisting of a charge-and-current preamplifier (PACI: Pre-Amplifier for Charge and Icurrent) [18] followed by sampling ADCs (digitizers). The preamplifiers had different gains and were located near the detectors, inside the vacuum chamber. The digitizers were in air, connected to the preamplifiers by differential cables. They were either 125 MS/s 12 bit cards for CsI [19] ("Florence cards" developed by INFN-Sezione di Firenze) or 100 MS/s 14 bit cards for silicons (developed by IPNO-Orsay, in the framework of the FAZIA collaboration).

The sampled waveforms delivered by the digitizers were then stored on disk and processed off-line in order to extract the relevant quantities for the analysis. The improvements in the identification capabilities presented in this paper are the result of several procedures which were developed during the R&D preparatory phase, see [4], [5], [6], [20], [21] and [22].

Identification techniques

The ∆E-E technique is based on the measurement of the energies deposited by the impinging particle in two detectors after punching through the first layer.

The PSA, on the other hand, requires the energy measurement of the reaction products stopped inside one detector and an additional parameter related to the shape of the charge signal or of the current signal [6]. In the present paper, this shape-related parameter is either the rise time of the charge signal for the "Energy vs Charge rise time" method, or the maximum of the current signal (called amplitude max in [22]) for the "Energy vs Current maximum" method. The algorithms and techniques for the digital signal processing employed in this analysis are described in [4] and [5]. Energy measurements as in [6].

During the first part of the experiment, particles impinged on the detectors through the high electric field side (front side injection) for both silicons (Si1 and Si2) of the telescope. Then for the second half of the experiment, we turned both the Si1 and Si2 by 180 o in order to compare the front and rear side injection configurations in the same conditions regarding detector characteristics, depletion voltages and electronics.

3 Rear or Front side particle injection? For the first silicon (Si1 of table 1), electric field (red lines) in V/m (red right scale) as a function of the depth for the two different mountings: front and rear side injection; particles enter from the left in both cases. The calculated energy losses [23], [24], [25] as a function of the range of the stopped particle in the semiconductor are also indicated in MeV/µm (black left scale) for alpha particles (for 2, 4 and 6 A MeV and they have been multiplied by a factor of 10 to scale the figure) and a 84 Kr nucleus (for 6, 12, 18 and 24 A MeV).

The time dependence of the current signal in semiconductors is ruled essentially by two phenomena: the charge carrier dissociation (the plasma erosion time) and the duration of the drift of the electrons (e -) and holes (h + ) towards the appropriate electrodes, which depends on electric field strength and carrier mobilities. The plasma erosion time is longer for higher densities of charge carriers along the incident ion ionization track, as it happens for higher atomic number Z. It depends also on the local electric field strength which varies with the injection side and penetration depth (see Fig. 1). For what concerns the charge carrier drift time, for a stopped particle in the first microns of the semiconductor (short track), in the rear side injection case, the electrons and holes experience first the low electric field. In particular the holes, which move towards opposite side (see Fig. 1), and have a mobility almost three times lower than the electrons, have to cross a major part of the thickness, thus increasing the collection time. In the front side injection case, for short tracks, particles stopped in the first few microns, the slower charge carriers (the holes) are collected by the close entrance side (front side) and their path lengths are shorter, thus reducing their collection time. At variance, for particles almost punching through the silicon layer (long tracks), the Bragg peak, responsible for the maximum density of charge carriers, stands near the exit side of the detector. So in the rear side injection case most of the charge carriers (those in the Bragg peak) "feel" a high electric field that strongly reduces the plasma erosion time. In the front side injection case the maximum density of charge carriers is located in the low electric field, thus plasma erosion time gets longer and signal rise times increase. This point is illustrated in Fig. 2 where charge signals for several nuclei stopped in the first silicon (Si1), at two different energies are shown in front (blue full lines) and rear (red dotted lines) side injection. The higher energy (larger signal amplitudes) corresponds to particles almost stopped at the end of the detector, whereas the low one refers to nuclei stopping just after the middle of the detector. We first observe that, at a given energy, the amplitude of the signal is independent of the entrance side. For particles almost punching through the first silicon junction, the signal shapes have a shorter rise time when the particles enter through the low electric field side (red dotted lines signal shapes of higher amplitude in each panel of Fig. 2. For nuclei stopped in the middle of the silicon layer (lower signals) the trend is opposite: rear side injection produces longer signal rise times. We recall that we are comparing results obtained with exactly the same detector and electronics. When the shape is used for discrimination, one may think that a larger signal shape variation would improve particle identification. Therefore, based on the shapes shown in Fig. 2, it intuitively appears that rear side injection will be, a priori, more favourable for PSA; it shows a larger variety of signal shapes and larger variations of the charge collection time.

Indeed in our previous experiment [6], as in other works [START_REF] Pausch | [END_REF], [7], [8], [10], [START_REF] Amorini | IEEE Nuclear Science Symposium Conference Record N[END_REF], [START_REF] Mutterer | [END_REF], ion discrimination was performed for particles entering the low electric field. We will now quantify the differences between the two mountings, by comparing the PSA performances (in terms of quality and energy thresholds) for rear and front side injection, using the two "Energy vs Charge rise time" and "Energy vs Current maximum" methods already used in the literature [START_REF] Pausch | [END_REF], [6], [10], [16]. However, we first want to explore the consequences of the silicon mounting on the ∆E-E identification technique (indeed it is mandatory to preserve first the excellent ∆E-E identification capabilities of the FAZIA detectors). So, in the next section, we will compare the ∆E-E performances for the two configurations. We expect, once the charge carriers are all collected and the signals properly filtered, both configurations to give the same quality of energy information: indeed, at variance with PSA, the ∆E-E method does not exploit the time development of the signals within the silicon.

4 Identification of fragments with ∆E(Si1)-E(Si2) method

In this section the performances of the ∆E-E method are compared in both configurations (rear and front side injection). In each silicon stage, the energy deposited is extracted from the total collected charge. We can see from Fig. 3 that, apart the slight difference due to the use of different targets (and thus slightly different cross sections for the production of elements), we obtain comparable identification maps in both cases, demonstrating the same quality of the energy measurement. A zoom at low energy for light nuclei, shown in Fig. 4, confirms the equally good isotopic identification at least for elements from H to Ti.

A more quantitative comparison is made in Fig. 5, where the previous maps have been linearized in order to get the Particle IDentification (PID) number distribution from neon to argon isotopes (similar results have been obtained for Z=1-9). Please note the logarithmic scale on the ordinate axis. The Gaussian isotopic distributions are completely separated at the level of their full width at half maximum height (FWHM), i.e. their FWHM never overlap for two adjacent isotopes. The quality of the isotopic separation is identical for front or rear injection, proving that for the ∆E-E method the energy measurement remains quantitatively the same, independent of the injection side. The only difference one may expect concerns the recombination of charge carriers along the track, contributing to Pulse Height Defect (PHD) phenomena that are quite important in the case of passivated ion-implanted silicon detectors [26].

Recombination leads to a loss in the collection of electrons and holes (and dN/dA a.u.

Z=18

Fig. 5. Linearization of the previous ∆E(Si1)-E(Si2) matrices, see Fig. 4, for both configurations, front (blue full line) and rear (red dashed line) side injection. For each atomic number Z, isotopic spectra have been normalized on the yields obtained for 21 Ne, 23 Na, 25 Mg, 27 Al, 29 Si, 32 P, 32 S, 35 Cl and 38 Ar respectively thus to a reduction of the measured energy) that is favoured by high charge densities and/or by a low electric field. Therefore this effect may depend on the injection side. The similarity of energy values measured for the two injection modes suggests that recombination is very small for the n-TD type detectors.

Comparison of Pulse Shape Analysis (PSA) in silicon for different injection modes

The PSA allows to identify particles stopped in one silicon detector from the information delivered by that detector only (Si1, 300 µm). Therefore, all results presented in this section have been obtained with a veto condition on the detector behind (Si2, 500 µm). The two methods of PSA used here, "Energy vs Charge rise time" and "Energy vs Current maximum", have been both successfully used by the FAZIA group and others, see [6], [16]. In order to have large dynamic ranges (≃2-3 GeV), low amplification values were used in the present experiment. A consequence of this choice was that no isotopic resolution with PSA could be observed, at variance with other FAZIA experiments using higher amplifications. Consequently the identification thresholds and Figure of Merit (FoM) [5], [6] calculated hereafter for PSA refer only to charge discrimination and are meant, in particular, to compare the two opposite injection faces.

I max (a.u.)

Rear injection

Fig. 6. Energy vs maximum amplitude of the current signal Imax for particles stopped in the first silicon (Si1). Particles punching through the detector have been removed.

PSA by energy vs maximum amplitude of the current signal

Figure 6 shows the "Energy vs Current maximum" correlation for both cases, front side injection (on the left) and rear side injection (on the right). For the rear side injection, we recover the identification map previously studied in the framework of the FAZIA collaboration, see [4], [5], [6]. The populated ridges allow to identify all atomic numbers from Z=1 to 36 (Kr beam). On the contrary, in the front side injection configuration no identification is visible. All elements merge together in a very compact cloud, corresponding to a strong correlation between the energy and the maximum current. Thus the maximum amplitude of the current signal is not a good PSA variable when the fragments enter through the high electric field side. We recall that only particles stopped in the first silicon of the telescope have been kept, the particles punching through have been removed.

PSA by energy vs rise time of the charge signal

Regarding the energy vs rise time of the charge signal shown in Fig. 7, we obtain in both cases identification maps similar to those presented in [5], [6] (for rear side injection) and [10] (for front side injection). As previously mentioned, we observe a very large increase of the maximum measured rise time (almost a factor of two, ≃200 ns instead of ≃100 ns) for the rear side injection configuration. Moreover, the shape of the correlation between energy and charge rise time is very different in the two cases. In the front side injection, the charge rise time continuously decreases with decreasing ion energy for all Z values. On the contrary, for rear side injection, starting from high kinetic energies, we observe a rise-and-fall trend of the rise time which, for slow ions, produces a ridge where all values merge together, whatever the particle is. In both cases a "no identification" zone is visible for each line at low energy, defining a Z dependent identification threshold. These thresholds will be determined more precisely in the following.

In Fig. 8 we present the same PSA matrices as in Fig. 7 but without removing the punching through particles, i.e. including particles which pass through Si1 and are stopped either in the second silicon or in the CsI. This presentation may be interesting because in some cases one might be tempted to apply the PSA in simply one silicon layer, without a veto on the second stage detector. This is the case of tagging systems mounted along the transport line of radioactive beams produced in fragmentation. The lack of a second layer qualitatively results in the spoiled correlations energy vs rise time of

Rear injection

Fig. 8. Energy vs rise time of the charge signal in the first silicon (Si1). Same as Fig. 7 but in this case punching through events have not been removed.

projectiles). From Fig. 8 one sees that in the rear side injection case punching through particles are located on the very left part of the identification matrix.

They are superposed on a rather narrow strip of high-energy stopped particles. Therefore, they introduce a threshold in the high energy zone of the PSA identification capability, but most of the PSA matrix is free from such a contamination. On the contrary, in the front side injection case, punching through particles lie just in the middle of the PSA matrix thus severely spoiling the identification. This constitutes a very strong limitation.

Particle identification thresholds for rear and front configuration

The PSA method based on "Energy vs Current maximum", does not allow to identify particles in the front side injection configuration. Therefore a comparison between front and rear side injection can be performed only on the basis of the second PSA technique, namely the "Energy vs Charge rise time" one. In this case, for the same stopped products, Fig. 7 shows very different identification maps. Even though the rear side injection method may seem more efficient, since it enlarges the ridges range, we need a more quantitative way to judge the identification quality for both configurations. Therefore we applied the same protocol as in [6]. A Figure of Merit (FoM) defined as

F oM = |P ID2 -P ID1| F W HM1 + F W HM2
was determined for adjacent Z as a function of the energy. Here FWHM1,2 are the full widths at half maximum of the Gaussian distributions of two adjacent elements of atomic number Z, Z+1. The corresponding PID1,2 stand for the peak values of the linearized ridges, as those obtained in case of the ∆E-E method in Fig. 5 for isotopic distributions. A value of FoM=0.7 was conventionally chosen in order to extract a low energy threshold above which we realize a good identification. A FoM value of 0.7 corresponds, for two peaks of equal intensity, to a peak-to-valley ratio of 2.0, see [5]. The same value FoM=0.7 leads to different energy thresholds of identification for different adjacent ejectiles. The method was applied to both matrices of Fig. 7. The identification thresholds are summarized in Fig. 9 in terms of energy (left) and energy per nucleon (right) using, for ions heavier than Z≥6, the mass parameterization A=2.08×Z+0.0029×Z 2 . As in [6], we have reported also the energy thresholds obtained with the standard ∆E-E technique (black symbols and line) inherently introduced by the first 300 µm thick silicon detector. One may clearly see a spectacular improvement on the identification energy threshold for the rear side injection technique (red symbols and line). For example, a gain of almost a factor of two in energy per nucleon is achieved for Z≃10-15. The thresholds obtained in the front side injection case are very similar to those published in [10] for other silicon junctions, indicating that they do not depend on the detector. At variance the corresponding rise time depends on the detector thickness and the applied voltage (electric field). Expressed as a function of the range in silicon, see Fig. 10, the threshold differences are even more impressive. In the front side injection case the mini- mum range for identification exhibits a concave shape with thresholds always larger than 150 µm, whereas in the rear side injection case the minimum range presents a continuous increase, from 30 to 160 µm. As already mentioned in [6] these thresholds do not depend on the detector thickness but it rather seems that they are a characteristic of reverse polarized silicon junctions [27]. Comparing the present results to those of [6], it also appears that the Z-identification thresholds increase with the bulk resistivity inhomogeneity.

Conclusions

In this experiment we have carefully compared the results of the digital "Pulse Shape Analysis"(PSA) technique for identifying stopped reaction products in a n-TD silicon detector for two different configurations: one in which the particles enter from the high electric field side (front side injection) and the other for incidence on the low electric field side (rear side injection). Two PSA techniques, previously successfully tested in the rear configuration, were examined. One, here called "Energy vs Current maximum", completely fails in the front side injection mode at variance with the rear side injection one. The other, "Energy vs Charge rise time", gives a very satisfactory result in both cases: full charge separation is obtained up to the maximum observed Z=36. However, in this latter case, the energy identification thresholds of atomic number Z are very different, with higher values for particles entering from the high electric field side. Among the two examined configurations, the present analysis indicates that the rear side injection mode has to be chosen for a better ion identification.

  Fig.1. For the first silicon (Si1 of table 1), electric field (red lines) in V/m (red right scale) as a function of the depth for the two different mountings: front and rear side injection; particles enter from the left in both cases. The calculated energy losses[23],[24],[25] as a function of the range of the stopped particle in the semiconductor are also indicated in MeV/µm (black left scale) for alpha particles (for 2, 4 and 6 A MeV and they have been multiplied by a factor of 10 to scale the figure) and a 84 Kr nucleus (for 6, 12, 18 and 24 A MeV).

Fig. 3 .Fig. 4 .

 34 Fig. 3. ∆E(Si1)-E(Si2) matrices for front (left panel) and rear (right panel) side injection. Particles punching through the second detector have been removed via a veto on the CsI signals.

Fig. 7 .

 7 Fig. 7. Energy vs rise time of the charge signal for particles stopped in the first silicon (Si1). Particles punching through the detector have been removed.

Fig. 8 (

 8 in our case most of the punching through ions are in the region of the Kr

Fig. 9 .

 9 Fig.9. Energy thresholds for Z identification with ∆E(300 µm)-E technique (black triangles) and with PSA technique (energy vs charge rise time: red points "rear side injection" and blue squares "front side injection") as a function of atomic number Z. The thresholds values are presented in terms of total energy (left) and energy per nucleon (right). Error bars take into account the statistic and the FoM method reproducibility to extract the energy thresholds.

Fig. 10 .

 10 Fig.10. Thresholds expressed in term of range in silicon material for Z identification with ∆E(300 µm)-E technique (black triangles) and with PSA technique (energy vs charge rise time: red points "rear side injection" and blue squares "front side injection") as a function of atomic number Z. Error bars take into account the statistic and the FoM method reproducibility to extract the energy thresholds.

Table 1

 1 Properties of the silicon detectors, CsI(Tl) and electronics used for ∆E-E and pulse shape measurements.

	Technical specifications	Si1	Si2	CsI
	Manufacturer	FBK	FBK	Amcrys
	Bulk type	n	n
	Thickness	300 µm	500 µm	10 cm+FBK photodiode
	Depletion voltage	132 V	290 V
	Applied voltage	140 V	260 V	30 V
	Resistivity unif.(FWHM)	≃4%	-
	Digitizer (bit/rate)	14 bit/100 MHz 14 bit/100 MHz	12 bit/125 MHz
	Exp. full Range	≥3 GeV	≃2 GeV	≃0.4 GeV

Acknowledgments The authors would like to thank the Superconducting Cyclotron staff of LNS, in particular D. Rifuggiato, for providing a very good quality beam. The support of the Machine shop team of LNS, G. Tobia and E. Scarlini of the Physics Department of Florence, S. Barbey from IPN (Orsay) and P. Desrue and Y. Merrer LPC (Caen), are gratefully acknowledged. Thanks are also due to N. Zorzi of FBK (Trento) for cooperative support in detector development. This work has been supported by the Agence Nationale de la Recherche under Grant number ANR-05-BLAN-0373-01. The research leading to these results has received funding from the European Union Seventh Framework Programme FP7 (2007-2013) under Grant Agreement n o 262010 -ENSAR and E.U. SPIRAL2 preparatory phase. We acknowledge support by the Foundation for Polish Science MPD program, co-financed by the European Union within the European Regional Development Fund.


Bernard Hanzon
email: b.hanzon@ucc.ie

Martine Olivi

Ralf Peeters

Subdiagonal pivot structures and associated canonical forms under state isometries

Many normalizations in various classes of linear dynamical state-space systems lead to system representations which are determined up to a state isometry. Here we present a new set of techniques to obtain (local) canonical forms under state isometries, using what we call sub-diagonal pivot structures. These techniques lead to a very flexible, straightforward algorithm to put any system into canonical form under state isometries. The parametrization of these canonical forms is discussed for a number of classes, including lossless systems and input-normal stable systems both in discrete time and in continuous time.

INTRODUCTION

Consider a linear time-invariant state-space system in discrete time with m inputs and p outputs:

x t+1 = Ax t + Bu t y t = Cx t + Du t t ∈ Z
or in continuous time ẋt = Ax t + Bu t y t = Cx t + Du t t ∈ R

x t ∈ R n for some nonnegative integer n (the state space dimension), u t ∈ R m and y t ∈ R p . The matrices A, B, C and D have real-valued entries and are of compatible sizes: n × n, n × m, p × n and p × m, respectively.

The corresponding transfer matrix of this system is given by G(z) = D+C(zI n -A) -1 B, which is a p×m matrix with rational functions as its entries. From realization theory it follows that, conversely, any p×m rational matrix function G(z) analytic at infinity can be written in the form

G(z) = D + C(zI n -A) -1 B
(1) where (A, B, C, D) is an appropriate quadruple of matrices and n a suitable state space dimension. Such a quadruple with the associated expression (1) is called a state-space realization of G(z).

The controllability matrix K and the observability matrix O associated with this system are defined as the blockpartitioned matrices

K = [B, AB, . . . , A n-1 B], O =     C CA . . . CA n-1     .
(2)

The system (or its input pair (A, B)) is called controllable if K has full row rank n and the system (or its output pair (C, A)) is called observable if O has full column rank n. where Q is an orthogonal matrix Q ′ Q = QQ ′ = I, and Q ′ denotes the transpose of Q.

Consider the following well-known characterizations of controllability:

(i) The pair (A, B) is controllable if and only if there exists T ∈ Gl n (R) such that T.K contains e 1 , e 2 , . . . , e n as columns, where e i denotes the ith standard basis vector e i = (0, 0, . . . , 0, 1, 0, . . . , 0) ′ with the number 1 in the i-th position. (ii) The pair (A, B) is controllable if and only if there exists T ∈ Gl n (R) such that (T.B, T.A.T -1 ) contains e 1 , e 2 , . . . , e n as columns in such a way that T.K contains e 1 , e 2 , . . . , e n as columns.

Consider the question whether there are, and if so what are analogous results for state isometries? Previously this has been investigated in [START_REF] Hanzon | Overlapping block-balanced canonical forms for various classes of linear systems[END_REF] and [START_REF] Peeters | Canonical lossless state-space systems: staircase forms and the schur algorithm[END_REF], where the structure of the controllability matrix plays a central role. Here we reconsider the question in a somewhat different way, where we stress the structure of the pair (A, B). This turns out to lead to much simpler answers.

Reasons for considering state isometries instead of state isomorphisms are:

v = (v(1), v(2), . . . , v(n)) ′ ∈ R n . The vector v is called a pivot vector with pivot at position k, 1 ≤ k ≤ n, or pivot-k vector for short, if v(k) > 0 and if v(j) with j > k are all zero. Definition 2.2. An n × r matrix M , r ≥ n is said to have a full pivot structure J = {j 1 , j 2 , . . . , j n } if for each k ∈ {1, 2, . . . , n} it holds that column j k of M is a pivot-k vector. Example 2.1. M =      * * * * * * + * + * * * * * 0 * 0 * * * + * 0 * 0 * + * 0 * 0 * 0 * 0 * 0 + 0 *     
where * denotes an arbitrary number and + denotes a (strictly) positive number, has a full pivot structure J = {7, 1, 5, 3, 6}. Definition 2.3. Consider the partitioned matrix [B|A] in R n× (m+n) . We say this has a sub-diagonal pivot structure if (i) [B|A] has a full pivot structure (ii) the prescribed pivot columns of A have the property that a column with pivot at position k has column number p k < k in A (hence column number

j k = m + p k < m + k in [B|A]), for each k ∈ {1, 2, . . . , n}. Example 2.2. [B|A] =    + * | * * * * 0 * | + * * * 0 + | 0 * * * 0 0 | 0 + * *    has a sub-diagonal pivot structure. Example 2.3. [B|A] =    + * | * * * * 0 * | + * * * 0 + | 0 * * * 0 0 | 0 * + *    has a sub-diagonal pivot structure. Example 2.4. [B|A] =    + * | * * * * 0 * | * * + * 0 + | * * 0 * 0 0 | * * 0 +   
does not have a sub-diagonal pivot structure. Theorem 1. If [B|A] has a sub-diagonal pivot structure then the pair (A, B) is controllable.

Proof. Recall Hautus (or PBH) criterion [START_REF] Hautus | Controllability and observability conditions of linear autonomous systems[END_REF]): the pair (A, B) is controllable if and only if ∄(w, λ) ∈ C n × C such that w * B = 0 and w * A = λw * , w = 0. Now suppose that w ∈ C n is such that w * B = 0 and there exists λ ∈ C such that w * A = λw * . Suppose further that w(i) = 0 and w(j) = 0 for all 1 ≤ j < i. Then, one of the following situations occurs (1) B contains the pivot-i vector:

w * B = 0 ⇔ w(i) = 0,
which yields a contradiction.

(2) A contains the pivot-i vector, say in column j, j < i:

w * A = λw * ⇔ w * Ae j = λw * e j = 0,
and since Ae j is the pivot-i vector, this again implies w(i) = 0 and we get a contradiction.

So by induction, w = 0. Theorem 2. For any given pivot structure which is not sub-diagonal, we can find an example [B|A] for which (A, B) is not controllable.

Proof. We distinguish two cases, depending on whether there is a pivot on the diagonal of the matrix A or not. Firstly, if there is a pivot on the diagonal position (i, i) of the matrix A then consider the case where all pivots are taken to be equal to 1 and all other entries of [B|A] are taken to be zero. Consider w = e i , the i-th standard basis vector. Clearly w ′ B = 0 because B does not contain the ith pivot vector. Furthermore w ′ A = e ′ i = w ′ because in its i-th row A has only one non-zero entry, namely the value 1 on the diagonal position (i, i) of A. Therefore according to the Hautus test, (A, B) is an uncontrollable pair.

Secondly, suppose there is a super-diagonal pivot. Let the super-diagonal pivot with the smallest row index be in position (i, j) in the matrix A, (hence in position (i, m + j) in the matrix [B|A]) where j > i. Let the j-th pivot be in column m + k in [B|A], where k ∈ {-m + 1, -m + 2, . . . , n}. Let [ B| Ã] denote the matrix that is obtained by taking all pivots equal to 1 and all other entries equal to zero. Now form the matrix

[B|A] = [ B| Ã] -E i,m+k + E i,m+i
where for any (r, s) with 1 ≤ r ≤ n, 1 ≤ s ≤ m + n, E r,s denotes the n × (m + n) matrix with all entries zero except for the entry (r, s), which is 1. Note that [B|A] has the same pivot structure as [ B| Ã], because in column m + k one has the j-th pivot and j > i, so the subtraction of E i,m+k does not change the pivot structure; and in column m + i there is no pivot on a super-diagonal position, because otherwise i would not be the smallest row-number with a super-diagonal pivot. Therefore column i either contains a sub-diagonal pivot or no prescribed pivot at all. In both cases addition of E i,m+i does not alter the pivot structure. Now note that w = e i + e j ∈ R n has the property that

w ′ [B|A] = [0|w ′ ]. This follows from the facts that e ′ i B = 0, e ′ i à = e ′ j , so e ′ i [ B| Ã] = e ′ m+j and e ′ j [ B| Ã] = e ′ m+k , e ′ i (-E i,m+k + E i,m+i ) = -e ′ m+k + e ′ m+i . So w = 0 is a left eigenvector of A in the left kernel of B, hence (A, B)
is uncontrollable according to the Hautus criterion.

CANONICAL FORMS UNDER STATE ISOMETRIES

We give some examples in which a normalization of the system by a state isomorphism leaves us the freedom of applying an orthogonal transformation (i.e. a state isometry).

Example 3.1. Let K n = [B|AB| . . . |A n-1 B] ∈ R n×nm denote the finite controllability matrix of (A, B, C, D).
Suppose that it has rank n. We can find a nonsingular square matrix T such that (T

K n )(T K n ) ′ = I n . Then ( Ã, B, C, D) := (T AT -1 , T B, CT -1 , D
) is normalized in the sense that its finite controllability matrix has orthonormal rows. Now note that this property is kept even if we apply a state isometry to the system: (Q ÃQ ′ , Q B, CQ ′ , D) has finite controllability matrix with orthonormal rows as well, for any orthogonal matrix Q. Example 3.2. One can do the analogous normalization for the finite observability matrix of the system. Example 3.3. If we have a stable discrete time system, we can define the controllability Gramian of the system in the usual way:

W c = K ∞ K ′ ∞ = ∞ k=0 A k BB ′ (A ′ ) k .
This Gramian solves the Lyapunov-Stein equation

W c -AW c A ′ = BB ′ . We can normalize the system to be input-normal, i.e. such that W c = I. If (A, B, C, D) is input-normal so is (QAQ ′ ,
c + W c A ′ = -BB ′ . The system is continuous-time input-normal if W c = I n . Similarly, the system is continuous-time output-normal if W o = I n where A ′ W o + W o A = -C ′ C.
Again the freedom of applying an orthogonal transformation remains. Example 3.6. Another well-known normalization is obtained by requiring the Gramians to be equal: W c = W o (this can be done both in the discrete time case and in the continuous time case). This can be called Gramianbalanced (or balanced in the sense of [START_REF] Helmke | Optimization and Dynamical Systems[END_REF]). Again this normalization is possible for any stable system and leaves the freedom of applying a state isometry. Effect of a state isometry is the same on both Gramians:

W c → QW c Q ′ and, using (Q -1 ) ′ = Q, W o → QW o Q ′ . Example 3.7. (A, B, C, D) is called Frobenius-balanced (or Euclidean balanced) if min T ∈Gln(R) à 2 F + B 2 F + C 2 F + D 2 F is attained at T = I where à = T AT -1 , B = T B, C = CT -1 , D = D and M 2 F = Tr(M M ′
). This normalization is again welldefined up to a state isometry (cf the book of [START_REF] Helmke | Optimization and Dynamical Systems[END_REF]). Example 3.8. The so-called LQG-input normal, boundedreal input normal and positive-real input normal forms described in [START_REF] Hanzon | Overlapping block-balanced canonical forms for various classes of linear systems[END_REF] again leave the freedom of applying a state isometry.

In all these cases it makes sense to look for canonical forms under state isometry. We now want to show that the sub-diagonal pivot structures correspond to local canonical forms under state isometry, and these local canonical forms are covering all cases. Theorem 3. Suppose that [B|A] and [QB|QAQ ′ ] have the same sub-diagonal pivot structure J = {j 1 , j 2 , . . . , j n } and that Q is orthogonal. Then Q = I.

Proof. First note that the pivot-1 column must be a column of B (it cannot lie in A) and of QB, so j 1 ≤ m. Therefore the first column of Q is of the form λe 1 , λ > 0, hence must be e 1 ; therefore Q can be partitioned as

Q =     1 0 . . . 0 0 . . . Q n-1 0     .
Now assume (induction hypothesis) we have shown that Q is of the form

Q = I i 0 0 Q n-i , i ∈ {1, . . . , n -1}
If the pivot-(i + 1) vector is in B we can use the same argument to conclude that the induction hypothesis holds for i + 1. If the pivot-(i + 1) vector is in A then: partition A as

A = A 11 A 12 A 21 A 22 , A 11 : i × i; A 21 : (n -i) × i.
We get

QAQ ′ = I i 0 0 Q n-i A 11 A 12 A 21 A 22 I i 0 0 Q ′ n-i = A 11 A 12 Q ′ n-i Q n-i A 21 Q n-i A 22 Q ′ n-i
The (i + 1)st pivot lies in the first row of A 21 and in the first row of

Q n-i A 21 hence the first column of Q n-i is e 1 . So Q =     I i+1 0 . . . 0 0 . . . Q n-i-1 0    
The theorem now follows by induction.

To show that every controllable pair (A, B) can be mapped to a sub-diagonal pivot structure by an appropriate state isometry we will present a recursive algorithm "which cannot get stuck, if controllability holds" and which is rather insensitive to small perturbations. We will use the well-known fact that for any vector x = 0 in R n-i there exists an orthogonal matrix Q n-i such that

Q n-i x =     x 0 . . . 0     First Step: (A, B
) is controllable so B = 0. It follows that B has at least one non-zero column. Choose one such column. Heuristically one can take the largest column. Find Q n , Q n Q ′ n = I n which maps this column to a vector of the form λe 1 , λ > 0. Recursion Step: Suppose the pivot-k column j k has been determined for k = 1, 2, . . . , i, in a sub-diagonal manner

(j k < m + k). Partition B = B(i) B(n -i) where B(i) is i × m; and A = Â(i, i) * A(n -i, i) * where Â(i, i) is i × i. Note that 1. [ B(i)| Â(i, i)] has a full pivot structure. 2. [B(n -i)|A(n -i, i)] has at least i zero columns. 3. [B(n -i)|A(n -i, i)] = 0; this can be seen as follows: if [B(n -i)|A(n -i, i)] = 0 then B = B(i) 0 , A = Â(i, i) * 0 *
hence the controllability matrix can be partitioned as

K = B(i) Â(i, i) B(i) Â(i, i) 2 B(i) . . . 0 0 0 . . .
which has rank i < n. This contradicts the assumption that (A, B) is controllable, q.e.d. 4. Application of a state isometry of the form

Q = I i 0 0 Q n-i transforms [B(n -i)|A(n -i, i)] into [Q n-i B(n -i)|Q n-i A(n -i, i)]
Choose a non-zero column from [B(ni)|A(ni, i)].

Heuristically one can take the largest column. Find

Q n-i , Q n-i Q ′ n-i = I n-i
which maps this column to a vector of the form λe 1 ∈ R n-i , λ > 0. Then the corresponding vector (in R n ) in [B|A] is a pivot-(i + 1) vector.

Applying the recursion step for each i = 1, 2, . . . , n -1, we obtain a sub-diagonal pivot structure for [B|A]. So we can construct a state isometry that brings [B|A] into a sub-diagonal pivot structure. It is clear that the same subdiagonal pivot structure can still be used under sufficiently small perturbations (namely if the length of the chosen (non-zero) vectors in each step does not reduce to zero).

Remarks

• This also gives a controllability test:

If (A, B) is con- trollable a sub-diagonal pivot structure is obtained. If (A, B) is not controllable an orthogonal Q is found such that (QAQ ′ , QB) = Â(i, i) * 0 * , B(i) 0
where [ B(i)| Â(i, i)] has a sub-diagonal pivot structure. When applied to input-normal discrete-time stable systems, [B|A] has orthonormal rows. In that case the size of the pivots is less than or equal to 1, so one may have a good idea whether a pivot is "numerically zero." • Model reduction by truncation. If [B|A] is in subdiagonal pivot form and we truncate the last ni rows and columns, we again have a controllable pair in sub-diagonal pivot form. • The staircase pivot structures (see [START_REF] Peeters | Canonical lossless state-space systems: staircase forms and the schur algorithm[END_REF]) form a sub-class of sub-diagonal pivot structures. A staircase pivot structure has the nice property that it corresponds to a pivot structure in the controllability matrix. However the recursive algorithm is not guaranteed to work for the staircase pivot structures. See section 6.

SUBDIAGONAL PIVOT STRUCTURE AND SCHUR PARAMETRIZATION OF DISCRETE-TIME LOSSLESS SYSTEMS

A system is called discrete-time lossless if it is stable and its m × m transfer matrix G(z) is unitary for all complex z with |z| = 1. It is well-known (cf., e.g., Proposition 3.2 in [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF] and the references given there) that the realization matrix, i.e. the block-partitioned matrix

R = D C B A (3)
associated with (1), is a real balanced minimal realization matrix of a lossless system if and only if R is an orthogonal matrix and A is asymptotically stable. It then holds that

W c = W o = I n .
For further background on lossless systems, see e.g. [START_REF] Genin | On Σ-lossless transfer functions and related questions[END_REF]. Note that there is a one-to-one correspondence between McMillan degree n lossless functions G(z), up to a left orthogonal factor, and controllable input-normal pairs [B|A] up to a state isometry.

In [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF] an atlas of overlapping balanced (local) canonical forms for lossless discrete-time systems of order n is presented: balanced state space realizations are constructed recursively, in line with the tangential Schur algorithm. Each such balanced canonical form is characterized by a fixed sequence of n numbers w k , |w k | < 1, k = 1, . . . , n, called the interpolation points, and a fixed sequence of n unit vectors u k ∈ R m , u k = 1, k = 1, . . . , n, called the direction vectors (which are not to be confused with the input signal applied to a system).

Here we will consider the case w k = 0, k = 1, . . . , n hence each balanced canonical form that we consider is determined by the choice of direction vectors. Each such balanced canonical form is then parameterized by an m × m orthogonal matrix D 0 and a sequence of n vectors v k , v k < 1, k = 1, . . . , n which are called the Schur parameter vectors.

In fact the orthogonal realization matrix can be written as a product of matrices of size

(m + n) × (m + n): R = Γ n Γ n-1 • • • Γ 1 Γ 0 ∆ T 1 ∆ T 2 • • • ∆ T n , (4
) where for k = 1, . . . , n:

Γ k = I n-k 0 0 0 V k 0 0 0 I k-1 , ∆ k = I n-k 0 0 0 U k 0 0 0 I k-1
with an (m+1)×(m+1) orthogonal matrix block V k given by

V k = v k I m -(1 -1 -v k 2 v k v T k v k 2 1 -v k 2 -v T k ,
and an (m + 1) × (m + 1) orthogonal matrix block U k given by

U k = u k I m -u k u T k 0 u T k and furthermore Γ 0 = I n 0 0 D 0 .
Note that here we consider the real case with real direction vectors and real Schur parameter vectors. Note further that Γ 0 , . . . , Γ n and ∆ 1 , . . . , ∆ n are all orthogonal matrices. It is important to note and not too difficult to see that the orthogonal matrix product

Γ = Γ n Γ n-1 • • • Γ 1 Γ 0
(5) in fact forms a positive m-upper Hessenberg matrix, i.e. an (m + n) × (m + n) matrix of which the m-th sub-diagonal only has positive entries and of which the last n -1 subdiagonals are all zero. It also follows almost directly that if the direction vectors u 1 , . . . , u n are taken to be standard basis vectors, then the matrix product

∆ T = ∆ T 1 ∆ T 2 • • • ∆ T n (6)
yields a permutation matrix. Hence in that case the balanced realization matrix R is obtained as a column permutation of an orthogonal positive m-upper Hessenberg matrix. More precisely, we have: Theorem 4. Let R be given by ( 4) and let the direction vectors be standard basis vectors: for k = 1, . . . , n, u n-k+1 = e i(k) , i(k) ∈ {1, . . . , m}. Then the sub-matrix [B|A] has a sub-diagonal pivot structure. This pivot structure is completely determined by the sequence i(1), i(2), . . . , i(n) as follows: Let k ∈ {1, 2, . . . , n}.

(i) If for all j < k, i(j) = i(k), then j k = i(k) and pivot-k is in B; (ii) otherwise let l = sup{j| j < k, i(j) = i(k)}, then j k = l + m and pivot-k is in the first k -1 columns of A.
Proof. Indeed, R = Γ∆ T in which Γ, given by ( 5), is a positive m-upper Hessenberg matrix and ∆ T the permutation matrix (6). The sub-matrix [B|A] is obtained from H = [0 I n ] Γ, which possesses the full pivot structure J = {1, 2, . . . , n}, by applying the permutation matrix ∆ T . For any j ∈ {1, 2, . . . , n}, post-multiplication by ∆ T n-j+1

only acts on the columns j to j + m as follows (where (i),(ii),(iii) apply simultaneously): (i) column j is moved into column position j + i(j) -1 (ii) column j + i(j) is moved into column position j + m (iii) columns j + 1, j + 2, . . . j + i(j) -1 and j + i(j) + 1, j + i(j) + 2, j + m are moved one position to the left, into column positions j, j +1, . . . j +i(j)-2 and j +i(j), j + i(j) + 1, j + m -1. Now consider k ∈ {1, 2, . . . , n}. The k-th column of H with a pivot at position k is only affected by the last k permutation matrices ∆ T n-k+1 , ∆ T n-k+2 , . . . , ∆ T n . In particular, by considering the special case k = 1, it follows that only the matrix ∆ T n acts on the first column of H (with pivot-1) which is moved into column i(1) of B. More generally, column k of H ends up in the first m + k -1 columns of [B|A], which ensures that [B|A] does have a sub-diagonal pivot structure.

The progression of column k of H with pivot at position k can be monitored as follows: For j = k : Post-multiplication by ∆ T n-k+1 moves column k into column k -1 + i(k). For each j = k -1, k -2, . . . as long as i(j) = i(k), postmultiplication by ∆ T n-j+1 moves the pivot-k vector one position to the left from column position j + i(k) into column position j -1 + i(k). Now, one of the following possibilities occurs: (i) for all 1 ≤ j < k, i(j) = i(k). Then column k of H is ultimately moved into column i(k) of B which inherits a pivot in position k.

(ii) there exists l, l = sup{j|k > j ≥ 1; i(j) = i(k)}, then eventually post-multiplication by ∆ T n-(l+1)+1 moves the pivot-k vector from column position l + i(l) into column position l+m where it remains and the pivot-k vector ends up in column l of A.

Conversely, we have Theorem 5. Let J = (j 1 , j 2 , . . . , j n ) be a given subdiagonal pivot structure for n × (m + n) matrices. For each k = 1, 2, . . . , n, choose the direction vector u n+1-k = e i(k) in the Schur algorithm by induction : i(1) := j 1 , and for k = 2, . . . , n

(i) if j k ≤ m, then i(k) := j k , (ii) if m < j k ≤ m + n define p k := j k -m < k (sub- diagonal structure); then i(k) := i(p k ).
For any choice of the Schur parameter vectors v 1 , v 2 , . . . , v n (all of length < 1) and for any choice of the orthogonal matrix D 0 , consider the (m + n) × (m + n) orthogonal realization matrix R given by (4). Then, the sub-matrix [B|A] possesses the sub-diagonal pivot structure J.

Remark. If a given lossless system allows a sub-diagonal pivot structure J then the previous result shows how to choose the corresponding direction vectors. Once these are known one can construct the corresponding Schur parameter vectors of the lossless system (cf [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF]). This implies that we have an explicit parametrization of the local canonical form of lossless systems that allow a given sub-diagonal pivot structure.

SUB-DIAGONAL PIVOT STRUCTURE AND CONTINUOUS-TIME LOSSLESS SYSTEMS

In the remark at the end of section 4 in [START_REF] Hanzon | Overlapping block-balanced canonical forms for various classes of linear systems[END_REF] on stable all-pass systems (i.e. on lossless systems) it was concluded that one obtains an atlas of local balanced canonical forms for lossless systems using the particular pivot structures presented in that paper. These pivot structures are staircase pivot structures hence are sub-diagonal. The atlas as well as the parametrizations presented there for each of the charts can be generalized straightforwardly to the case where one allows all possible sub-diagonal pivot structures. The parametrization of a local balanced canonical form associated to a sub-diagonal pivot structure J = (j 1 , j 2 , . . . , j n ) is obtained as follows:

Firstly let [B| Ã] be an n × (m + n) (real) matrix with sub-diagonal pivot structure J and with skew-symmetric matrix Ã, i.e. Ã + ÃT = 0. The set of all such matrices can be parameterized by R n + ×R (m-1)n where the positive numbers are describing the pivots and the other numbers are describing the remaining free entries in B and in the strictly lower triangular part of the skew symmetric matrix Ã. Secondly let V be the unique upper triangular matrix satisfying V + V T = -BB T . Thirdly choose D to be an arbitrary m × m orthogonal matrix and let C = -DB T . Finally let A = Ã + V. Then (A, B, C, D) is the balanced realization with the sub-diagonal pivot structure J, of a lossless system. Note that A + A T = -BB T = -C T C showing that the state-space system is balanced. In this way one obtains a chart for each sub-diagonal pivot structure J; varying over all sub-diagonal pivot structures one obtains an atlas covering the manifold of McMillan degree n continuous-time lossless systems.

ON THE RELATION OF SUB-DIAGONAL PIVOT STRUCTURES WITH STAIRCASE FORMS

As mentioned previously, the atlas studied in [START_REF] Peeters | Canonical lossless state-space systems: staircase forms and the schur algorithm[END_REF] provides overlapping balanced canonical forms with a particular sub-diagonal pivot structure : the so-called staircase forms. In these forms, the matrix [B|A] has a particular sub-diagonal pivot structure: if the columns of B contain p B pivots, the remaining p A = np B pivots have to be located in the first p A columns of A with increasing pivot positions. This distribution of pivots implies that the associated controllability matrix contains a pivot at position k for k = 1, 2, . . . , n so the controllability matrix has a full pivot structure.

The following representation of a sub-diagonal pivot structure J = (j 1 , j 2 , . . . , j n ) allows to clarify the connection. An m × n matrix Y = (y i,j ) can be associated with the sub-diagonal pivot structure. Its entries are all zeros except n entries corresponding to the pivot positions 1, 2, . . . , n. This matrix is constructed by induction : Pivot-1 is in the j 1 -th column of B and we put Y j1,1 = 1 Assume that the pivot positions 1, 2, . . . , k -1 have been displayed, then (i) either pivot-k is in the j k -th column of B, and we put

Y j k ,1 = k (ii) or pivot-k is in the p k -th column of A. Since [B|A]
has a sub-diagonal structure, then p k ≤ k -1 has already been displayed: Y i,j = p k for some i and j. Then we put Y i,j+1 = k.

In the case of a staircase form, this matrix has a special structure and specifies the pivot structure of the controllability matrix. In view of its relationship with nice selections, it is called a numbered Young diagram (see [Peeters et al., 2007, section 4]).

The recursive algorithm presented in the previous section provides a chart selection algorithm, that was lacking for the staircase forms. To be more precise: it is known that, if (A, B) is controllable, then, if one knows in advance which choices of pivots in B, can be made, one can, in order to obtain a staircase form, again apply the recursive algorithm; in that case one will have to choose the left-most possible pivot in A whenever a pivot is not chosen from B at any step in the recursive algorithm. The difficulty here is of course that usually one will not know in advance which choices of pivots from B lead to a staircase form.

If one makes a wrong choice for the pivots of B in the recursive algorithm, then at some stage in the recursive algorithm in which a pivot, say pivot i, from A needs to be chosen in column j of A, the relevant column vector (a i,j , a i+1,j , . . . , a n,j ) ′ will be equal to the zero vector in R n-i+1 . In that case one can still continue the recursive algorithm to obtain a sub-diagonal pivot structure.

If at each stage one takes either a pivot in B or one takes the left-most possible pivot in A, then the corresponding controllability matrix will have a pivot structure! Using that one can use the numbered Young diagram methods of [START_REF] Peeters | Canonical lossless state-space systems: staircase forms and the schur algorithm[END_REF] to determine a choice of pivots in B that correspond with a staircase form. This can then be used in a second round to obtain a staircase form using the recursive algorithm. This should still work under effects of (sufficiently small) round-off error or more general, when one decides to ignore (sufficiently) small perturbations.

CONCLUSION

Usage of sub-diagonal pivot structures and associated algorithms can lead to new flexible canonical forms, for which parametrizations are available in several cases, including discrete-time and continuous-time lossless systems. The results can also be used to construct controllability tests and procedures for model reduction by truncation. As only orthogonal transformations are used, procedures will be numerically stable. Generalizations to the complex-valued case are straightforward.

  Minimality holds iff both controllability and observability hold, which holds iff the McMillan degree of G(z) is equal to n. Two minimal realizations (A, B, C, D) and ( Ã, B, C, D) of a given function G(z) are always similar: there exists a unique invertible matrix T such that ( Ã, B, C, D) = (T AT -1 , T B, CT -1 , D). We define a state isomorphism as a map : (A, B, C, D) → (T AT -1 , T B, CT -1 , D), T ∈ Gl n (R) and a state isometry as a map (A, B, C, D) → (QAQ ′ , QB, CQ ′ , D),

  QB, CQ ′ , D) for any orthogonal Q. Example 3.4. One can do the analogous normalization using the observability Gramian W o . It satisfies W o -A ′ W o A = C ′ C. The resulting system is output-normal if W o = I. Again there is the freedom of applying a state isometry without losing this property. Example 3.5. Assume the eigenvalues of A are in the open left half-plane. Then we can define the continuous time controllability Gramian W c . It satisfies the Lyapunov equation AW

•

  Application of orthogonal group elements is in general numerically more stable than application of elements of the general linear group Gl n (R).• State isometries appear in the theory of input-normal forms, output-normal forms, "equal Gramians" balanced forms (where the Gramians do not have to be diagonal, just equal), balanced forms for lossless

	systems, and other normal forms for various classes
	of linear systems (some examples are given in section
	3).
	2. PIVOT STRUCTURES IN INPUT PAIRS
	Definition 2.1. Let n be a positive integer. Consider a
	vector