Identification of single-input-single-output quantum linear systems

The purpose of this paper is to investigate system identification for single-input-single-output general (active or passive) quantum linear systems. For a given input we address the following questions: (1) Which parameters can be identified by measuring the output? (2) How can we construct a system realization from sufficient input-output data? We show that for time-dependent inputs, the systems which cannot be distinguished are related by symplectic transformations acting on the space of system modes. This complements a previous result for passive linear systems. In the regime of stationary quantum noise input, the output is completely determined by the power spectrum. We define the notion of global minimality for a given power spectrum, and characterize globally minimal systems as those with a fully mixed stationary state. We show that in the case of systems with a cascade realization, the power spectrum completely fixes the transfer function, so the system can be identified up to a symplectic transformation. We give a method for constructing a globally minimal subsystem direct from the power spectrum. Restricting to passive systems the analysis simplifies so that identifiability may be completely understood from the eigenvalues of a particular system matrix.

System identification theory [32][33][34][35][36] lies at the interface between control theory and statistical inference, and deals with the estimation of unknown parameters of dynamical systems and processes from input-output data. The integration of control and identification techniques plays an important role, e.g., in adaptive control [37]. The identification of linear systems is by now a well developed subject in classical systems theory [32][33][34][38][39][40][41][42][43][44][45], but has not been fully explored in the quantum domain [1]. This paper deals with the problem of identifying unknown dynamical parameters of quantum linear systems (QLSs). A QLS is a continuous variables open system with modes a = (a 1 , . . . , a n ) T , which has a quadratic * Electronic address: pmxml2@nottingham.ac.uk is the infinitesimal annihilation operator at time t. The joint dynamics is completely char-arXiv:1608.01227v3 [quant-ph] 22 Dec 2017 acterized by the triple (S, C, Ω) consisting of a 2m × 2m scattering matrix S, a 2m × 2n system-input coupling matrix C, and a 2n × 2n Hamiltonian matrix Ω. Since each system or channel mode has two coordinates corresponding to creation and annihilation operators, all matrices have a 2 × 2 block structure, and it is convenient to use the "doubled-up" conventions introduced in [17], as detailed in Sec. II. The data (S, C, Ω) fix the joint unitary dynamics U(t) obtained as a solution of a quantum stochastic differential equation [46]; due to the quadratic interactions, the evolved modes a(t) := U(t) † aU(t) and output fields B out (t) := U(t) † B(t)U(t) are linear transformations of the original degrees of freedom.
In a nutshell, system identification deals with the estimation of dynamical parameters of input-output systems from data obtained by performing measurements on the output fields. We distinguish two contrasting approaches to the identification of linear systems, which we illustrate in Fig. 1. In the first approach, one probes the system with a known time-dependent input signal (e.g., coherent state), then uses the output measurement data to compute an estimator of the unknown dynamical parameter. In the Laplace domain, the input and output fields are related by a linear transformation given by the 2m × 2m transfer function Ξ(s): whereb(s) is the vector of input creation and annihilation input noise operators. The transfer function Ξ(s) is a rational matrix valued function, which becomes a symplectic matrix in the "frequency domain" (i.e., for s = −iω ∈ iR), reflecting the fact that the unitary dynamics preserves the canonical commutation relations.
Similarly to the classical case, Eq. (1) means that the input-output data can be used to reconstruct the transfer function Ξ(s), while systems with the same transfer function cannot be distinguished. Therefore, the basic identifiability problem is to find the equivalence classes of systems with the same transfer function.
In [1] this problem was analyzed for the special class of passive quantum linear systems (PQLSs) and it was shown that minimal equivalent systems are related by n × n unitary transformations acting on the space of annihilation modes a. By definition a QLS is minimal if no lower dimensional system has the same transfer function, which in the passive case is equivalent to the system being either observable, controllable, or Hurwitz stable [1]. In Sec. III we answer the identifiability question for the case of general (not necessarily passive) QLSs; we show that the equivalence classes are determined by symplectic transformations acting on the doubled-up space of canonical variablesȃ. It is worth noting that while in the classical set-up equivalent linear systems are related by similarity transformations, in both quantum scenarios described above the transformations are more restrictive due to the unitary nature of the dynamics.
In the second approach, the input fields are prepared in a stationary in time, pure Gaussian state with inde-pendent increments (squeezed vacuum noise), which is completely characterised by the covariance matrix V = V (N, M ) and the associated quantum Ito rule [17] If the system is minimal and Hurwitz stable, the dynamics exhibits an initial transience period after which it reaches stationarity and the output is in a stationary Gaussian state, whose covariance in the frequency domain is given by the power spectrum Since the power spectrum depends quadratically on the transfer function, the parameters which are identifiable in the stationary scenario will also be identifiable in the time-dependent one. Our goal is to understand to what extent the converse is also true. First, we note that for a given minimal system there may exist lower dimensional systems with the same power spectrum. To understand this, consider the system's stationary state and note that it can be uniquely written as a tensor product between a pure and a mixed Gaussian state (cf. the symplectic decomposition). In Theorem 2 we show that restricting the system to the mixed component leaves the power spectrum unchanged. Furthermore, the pure component is passive, which ties in with previous results of [23]. Conversely, if the stationary state is fully mixed, there exists no smaller dimensional system with the same power spectrum. Such systems will be called globally minimal, and can be seen as the analog of minimal systems for the stationary setting.
One of the main results is Theorem 3 which shows that for "generic" globally minimal single-input-single-output (SISO) systems which admit a cascade representation, the power spectrum Ψ V (s) determines the transfer function Ξ(s) uniquely, and therefore the time-dependent and time-stationary identifiability problems are equivalent. It is interesting to note that this equivalence is a consequence of unitarity and purity of the input state, and does not hold for generic classical linear systems [38,41].
The paper is structured as follows. In Sec. II we review the setup of input-output QLSs, and their associated transfer function. We discuss in greater detail the two identifiability approaches mentioned above. In Sec. III we study the identifiability of QLSs in the timedependent input setting. In Theorem 1 we show that the equivalence classes of input-output systems with the same transfer function are given by symplectic transformations of the system's modes. We further show how a physical realization can be constructed from the system's transfer function. In Sec. IV we analyze the identifiability of QLSs in a stationary Gaussian noise input setting. We introduce the notion of global minimality for systems with minimal dimension for a given power spectrum, and show that a system is globally minimal if and only if it has a fully mixed stationary state, cf. Theorem 2. In Sec.
V we analyze the structure of the power spectrum identifiability classes, and show that the power spectrum determines the transfer function uniquely, for a large class of SISO systems, cf. Theorem 3. Finally, we show that using an additional input channel with an appropriately chosen entangled input ensures that the system is always globally minimal.

A. Preliminaries and notation
We use the following notations: "Tr" and "Det" denotes the trace and determinant of a matrix, respectively. For a matrix X = (X ij ) the symbols: X # = (X * ij ), X T = (X ji ), X † = (X * ji ) represent the complex conjugation, transpose, and adjoint matrix respectively, where "*" indicates complex conjugation. We also use the doubled-up notationX := X T , (X # ) T T and ∆(A, B) := A, B; B # , A # . For example, we may write is the set of all distinct eigenvalues of X. A similar notation is used for matrices of operators. We use "1" to represent the identity matrix or operator. δ jk is Kronecker δ and δ(t) is Dirac δ. The commutator is denoted by [·, ·].
If additionally, S is of the form S = ∆(S − , S + ) for some S − , S + ∈ R m×m then we say that it is symplectic. Such matrices form a group called the symplectic group [17,47].

II. QUANTUM LINEAR SYSTEMS
In this section we briefly review the QLS theory, highlighting along the way results that will be relevant for this paper. We refer to [26] for a more detailed discussion on the input-output formalism, and to the review papers [13,15,46,48] for the theory of linear systems.

A. Time-domain representation
A linear input-output quantum system is defined as a continuous variables (cv) system coupled to a Bosonic environment, such that their joint evolution is linear in all canonical variables. The system is described by the column vector of annihilation operators, a := [a 1 , a 2 , . . . , a n ] T , representing the n cv modes. Together with their respective creation operators a # := [a # 1 , a # 2 , . . . , a # n ] T they satisfy the canonical commutation relations (CCR) a i , a * j = δ ij 1. We denote by H := L 2 (R n ) the Hilbert space of the system carrying the standard representation of the n modes. The environment is modelled by m bosonic fields, called input channels, whose fundamental variables are the fields T , where t ∈ R represents time. The fields satisfy the CCR Equivalently, this can be written as are the infinitesimal (white noise) annihilation operators formally defined as b i (t) := dB i (t)/dt [15]. The operators can be defined in a standard fashion on the Fock space F = F(L 2 (R) ⊗ C m ) [5]. For most of the paper we consider the scenario where the input is prepared in a pure, stationary in time, meanzero, Gaussian state with independent increments characterized by the covariance matrix where the brackets denote a quantum expectation. Note that N = N † , M = M T , and V ≥ 0, which ensures that the state does not violate the uncertainty principle. The state's purity can be characterized in terms of the symplectic eigenvalues of V , as will be discussed in Sec. IV. In particular, N = M = 0 corresponds to the vacuum state, while pure squeezed states for single-input-single-output (SISO) systems (i.e., m = 1) satisfy |M | 2 = N (N + 1). More generally, we consider a nonstationary scenario where the input state has time-dependent mean B(t) , e.g., a coherent state with time-dependent amplitude. For more details on Gaussian states see [49,50]. The dynamics of a general input-output system is determined by the system's Hamiltonian and its coupling to the environment. In the Markov approximation, the joint unitary evolution of system and environment is described by the (interaction picture) unitary U(t) on the joint space H ⊗ F, which is the solution of the quantum stochastic differential equation [5,26,46,48,51] with initial condition U(0) = I. Here, H and L are system operators describing the system Hamiltonian and coupling to the fields; dB i (t), dB # i (t), are increments of fundamental quantum stochastic processes describing the creation and annihilation operators in the input channels.
For the special case of linear systems, the coupling and Hamiltonian operators are of the form for m × n matrices C − , C + and n × n matrices Ω − , Ω + satisfying Ω − = Ω † − and Ω + = Ω T + . As shown below, this ensures that all canonical variables evolve linearly in time. Indeed, let a(t) and B out (t) be the Heisenberg evolved system and output variables a(t) := U(t) † aU(t), B out (t) := U(t) † B(t)U(t). (5) By using the QSDE (4) and the Ito rules (3) one can obtain the following Ito-form quantum stochastic differential equation of the QLS in the doubled-up notation [17] whereȃ := (a T , a # T ) T , C := ∆ (C − , C + ), and A := ∆ (A − , A + ) = − 1 2 C C − iJ n Ω with Ω = ∆ (Ω − , Ω + ) and It is important to note that not all choices of A and C may be physically realizable as open quantum systems [16]. A special case of linear systems is that of passive quantum linear systems (PQLSs) for which C + = 0 and Ω + = 0, whose system identification theory was studied in [1]. We will return to this important class along the way. This type of system often arises in applications, and includes optical cavities and beam splitters.

B. Controllability and observability
By taking the expectation with respect to the initial joint system state of Eqs. (6) we obtain the following classical linear system Definition 2. The quantum linear system (6) is said to be Hurwitz stable (respectively controllable, observable) if the corresponding classical system (8) is Hurwitz stable (respectively controllable, observable).
In general, for a quantum linear system observability and controllability are equivalent [21]. A system possessing one (and hence both) of these properties is called minimal. Checking minimality comes down to verifying that the rank of the following observability matrix is 2n: where Ω = ∆(Ω − , Ω + ). In the case of passive systems Hurwitz stability is further equivalent to minimality of the system [1]. However for active systems, although the statement [Hurwitz =⇒ minimal] is true [23], the converse statement ([minimal =⇒ Hurwitz]) is not necessarily so. We see this by means of a counterexample. Example 1. Consider a general one-mode SISO QLS, which is parametrizsed by Ω = ∆(ω − , ω + ) and C = ∆ (c − , c + ). The system is Hurwitz stable (i.e. the eigenvalues of A have a strictly negative real part) if and only if (1) |c − | > |c + | and |ω − | ≥ |ω + |, or A system is nonminimal if and only if the following matrix has rank less than 2: Clearly it is possible for a system to be {minimal}∩{Hurwitz} or {non-minimal}∩{non-Hurwitz}.
Further, for a counterexample to the statement: [minimal =⇒ Hurwitz] consider for example In light of the previous example, we make the physical assumption that all systems considered throughout this paper are Hurwitz (hence minimal).

C. Frequency-domain representation
For linear systems it is often useful to switch from the time domain dynamics described above, to the frequency domain picture. Recall that the Laplace transform of a generic process x(t) is defined by where s ∈ C. In the Laplace domain the input and output fields are related as follows [9]: where Ξ(s) is the transfer function matrix of the system In particular, the frequency domain input-output re- and similarly for the output modes [53]. As a consequence, the transfer matrix Ξ(−iω) is symplectic for all frequencies ω [17].
More generally one may allow for static scattering (implemented by passive optical components such as beamsplitters) or static squeezing processes to act on the interacting field before interacting with the system. The corresponding transfer function is obtained by multiplying the transfer function (12) with the scattering or squeezing symplectic matrix S on the right [17].
In the case of passive systems, Ξ + (s) ≡ 0 and so the doubled-up notation is no longer necessary; the inputoutput relation becomes [1,9] b out (s) = Ξ(s)b(s), (13) where the transfer function is given by which is unitary for all s = −iω ∈ iR. In the case of passive systems we write the triple determining the evolution as (S, C − , Ω − ), where the scattering matrix S is unitary. Finally, we note that while the transfer function is uniquely determined by the triple (S, C, Ω), the converse statement is not true, as discussed in detail in the next section.

A. Identifiability classes
We now consider the following general question: which dynamical parameters of a QLS can be identified by observing the output fields for appropriately chosen input states? This is the quantum analog of the classical system identification problem addressed in [39][40][41]. The inputoutput relation (11) shows that the experimenter can at most identify the transfer function Ξ(s) of the system. Systems which have the same transfer function are called equivalent and belong to the same equivalence class.
Before answering this question for general QLSs we discuss the case of passive QLSs considered in [1]. The transfer function in Eq. (13) can be identified by sending a coherent input signal of a given frequency ω and known amplitude α(ω), and measuring the output state, which is a coherent state of the same frequency and amplitude Ξ(−iω)α(ω).
In the case of passive systems it is known that two minimal systems with parameters (Ω, C, S) and (Ω , C , S ) are equivalent if and only if their parameters are related by a unitary transformation, i.e. C = CT and Ω = T ΩT † for some n×n unitary matrix T , and S = S . The first part of this result was shown in [1]; the fact that the scattering matrices must be equal follows by choosing s = −iω and taking the limit ω → ∞ in Eq. (14). Physically, this means that at frequencies far from the internal frequencies of the system, the input-output is dominated by the scattering or squeezing between the input fields. Our first main result is to extend this result to general (active) linear systems.
Proof. Firstly, using the same argument as above, the scattering or squeezing matrices S and S must be equal. It is known [32] that two minimal classical linear systems for input u(t), output y(t), and system state x(t) have the same transfer function if and only if Note that at this stage T is not assumed to be symplectic. The second and third conditions imply C = C T T , which further implies that [T T, C C] = 0. Now by earlier definitions A = − 1 2 C C − iJ n Ω, so that the second and third conditions applied to the first condition imply that J n Ω = T J n ΩT −1 . Next, using this and the observation (J n Ω) = J n Ω it follows that [T T, J n Ω] = 0. Now, C (J n Ω) k = C T T (J n Ω) k = C (J n Ω) k T T which means that the minimality matrix O satisfies O = OT T . Because the system is minimal O must be full rank, hence T T = 1.
Finally, it remains to show that the matrix T generating the equivalence class is of the form To see this, observe that CA k , C A k must be of the of this doubled up form for k ∈ {0, 1, 2, . . .}. Writing CA k , C A k , and T as T3 T4 , and using the above result, C A k = CA k T , it follows that and so using the fact that O is full rank gives the required result.
Therefore, without any additional information, we can at most identify the equivalence class of systems related by a symplectic transformation (on the system). Note that the above transformation of the system matrices is equivalent to a change of co-ordinatesȃ → T ȃ in Eq. (6).

B. Identification method
Suppose that we have constructed the transfer function from the input-output data, using for instance one of the techniques of [32] and [54].
Here we a outline a method to construct a system realization directly from the transfer function, for a general SISO quantum linear system. The realization is obtained indirectly by first finding a non-physical realization and then constructing a physical one from this by applying a criterion developed in [21]. The construction follows similar lines to the method described in [1] for passive systems.
Let (A 0 , B 0 , C 0 ) be a triple of doubled-up matrices which constitute a minimal realization of Ξ(s), i.e., For example, in Appendix A such a realization is found for an n-mode minimal SISO system, with matrices (A, C), possessing 2n distinct poles each with a non-zero imaginary part. Any other realization of the transfer function can be generated via a similarity transformation The problem here is that in general these matrices may not describe a genuine quantum system in the sense that from a given A, B, C one cannot reconstruct the pair (Ω, C). Our goal is to find a special transformation T mapping (A 0 , B 0 , C 0 ) to a triple (A, B, C) that does represent a genuine quantum system. Such triples are characterized by the following physical realizability conditions [21] A + A + C C = 0 and B = −C .
Therefore, substituting (17) into the left equation of (18) one finds where the matrices J here are of appropriate dimensions. Next, because the system is assumed to be stable it follows from [ We now need to use a result from [19], which is a sort of singular value decomposition for symplectic matrices. We state the result in a slightly different way here. Lemma 1. Let N 2n×2n be a complex, invertible, doubledup matrix and let N = N N .
(1) Assume that all eigenvalues of N are semisimple [55]. Then there exists a symplec- Here is one of the Pauli matrices and 1 2 is the identity.
(2) There exists another symplectic matrix V such that The coefficients α i and β i are determined from µ i and ν i via The lemma can be extended beyond the semisimple assumption, but since the latter holds for generic matrices [19], it suffices for our purposes.
We can therefore use Lemma 1 together with Eq. (20) in order to write the "physical" T as T = VT W , where W andT can be computed as in the lemma above, and V is a symplectic matrix. However, since the QLS equivalence classes are characterized by symplectic transformation, this means that T 0 =T W transforms (A 0 , B 0 , C 0 ) to the matrices of a quantum systems satisfying the realizability conditions. Finally, we can solve to find the set of physical parameters (Ω, C), which are given in terms of (A 0 , B 0 , C 0 ), as Remark 1. Note that, by assumption, Ξ(s) is the transfer function of a QLS. Since the original triple (A 0 , B 0 , C 0 ) is minimal, this implies that there exists a nonsingular T satisfying (20), so the right side of (20) is nonsingular, which eventually leads to a nonsingular transformation T computed using Lemma 1.
Remark 2. The proof also holds for multiple-inputmultiple-output (MIMO) systems provided that one can find a minimal doubled-up (non-physical) realization beforehand.

C. Cascade realization of QLS
Recently, a synthesis result has been established showing that the transfer function of a "generic" QLS has a pure cascade realization [18]. Translated to our setting, this means that given a n-mode QLS (C, Ω), one can construct an equivalent system (i.e., with the same transfer function) which is a series product of single mode systems. The result holds for a large class of systems characterized by the fact that the matrix A admits a certain symplectic Schur decomposition, which holds for a dense, open subset of the relevant set of matrices.
Assuming that such a cascade is possible, the transfer function is an n-mode product of single mode transfer functions, which are given by Further, we can stipulate that the coupling to the field is of the form C = ∆(C − , 0), with each element of C − being real and positive. Indeed, since the system is assumed to be stable, there exists a local symplectic transformation on each mode so that coupling is purely passive. The point of this requirement is that it fixes all the parameters, so that under these restrictions each equivalence class from Sec. III contains exactly one element. Note that the Hamiltonian may still have both active and passive parts. Therefore, each one mode system in the series product is characterized by three parameters, c i , Ω i− ∈ R with c i = 0, and Ω i+ ∈ C. If Ω i+ = 0 then the mode is passive. Actually, it is more convenient for us here to reparametrize the coefficients so that . Therefore, from the properties of the individual Ξ i± (s), one finds that Ξ − (s) and Ξ + (s) can be written as with γ, γ i , λ i ∈ C, x i ∈ R, and y i either real or imaginary, while j is some number between 1 and n−1. In particular, the poles are either in real pairs or in complex-conjugate pairs. Furthermore, there is a possibility that some of the poles and zeros may cancel in (21) and (22), and as a result some of these poles and zeros could be fictitious (see proof of Theorem 3 later where this becomes important).
For passive systems such a cascade realization is always possible [14,21] and each single mode system is passive. We show how this may be done in the following example.
Example 2. Consider a SISO PQLS (C, Ω) and let z 1 , z 2 , . . . , z m be the eigenvalues of A = −iΩ − 1 2 C † C. Then the transfer function is given by Now, comparing each term in the product with the transfer function of a SISO system of one mode, i.e., it is clear that each represents the transfer function of a bona-fide PQLS with Hamiltonian and coupling parameters given by Ω i = −Im(z i ) and 1/2|c i | 2 = −Re(z i ). This realization of the transfer function is a cascade of optical cavities. Furthermore, we note that the order of the elements in the series product is irrelevant; in fact a differing order can be achieved by a change of basis on the system space (see Sec. III).
In actual fact this result enables us to find a system realization directly from the transfer function, thus offering a parallel strategy to the realization method in Sec. III B for passive systems. Note that a similar brute-force approach for finding a cascade realization of a general SISO system is also possible. However, the active case is more involved than the passive case, as the transfer function is characterized by two quantities, Ξ − (s) and Ξ + (s), rather than just one. For this reason and also that Sec. III B indeed already offers a viable realization anyway, we do not discuss the result here.

IV. POWER SPECTRUM SYSTEM IDENTIFICATION
Until now we addressed the system identification problem from a time-dependent input perspective. We are now going to change viewpoint and consider a setting where the input fields are stationary (quantum noise) but may have a non-trivial covariance matrix (squeezing). In this case the characterization of the equivalence classes boils down to finding which systems have the same power spectrum, a problem which is well understood in the classical setting [41] but has not been addressed in the quantum domain.
The input state is "squeezed quantum noise", i.e., a zero-mean, pure Gaussian state with time-independent increments, which is completely characterized by its covariance matrix V = V (M, N ) cf. Eq. (3). In the frequency domain the state can be seen as a continuous tensor product over frequency modes of squeezed states with covariance V (M, N ). Since we deal with a linear system, the input-output map consists of applying a (frequency dependent) unitary Bogoliubov transformation whose linear symplectic action on the frequency modes is given by the transfer function Consequently, the output state is a Gaussian state consisting of independent frequency modes with covariance matrix where Ψ V (−iω) is the restriction to the imaginary axis of the power spectral density (or power spectrum) defined in the Laplace domain by Our goal is to find which system parameters are identifiable in the stationary regime where the quantum input has a given covariance matrix V . Since in this case the output is uniquely defined by its power spectrum Ψ V (s) this reduces to identifying the equivalence class of systems with a given power spectrum. Moreover, since the power spectrum depends on the system parameters via the transfer function, it is clear that one can identify "at most as much as" in the time-dependent setting discussed in Sec. III. In other words the corresponding equivalence classes are at least as large as those described by symplectic transformations (15).
In the analogous classical problem, the power spectrum can also be computed from the output correlations. The spectral factorization problem [42] is tasked with finding a transfer function from the power spectrum. There are known algorithms [42,44] to do this. From the latter, one then finds a system realization (i.e. matrices governing the system dynamics) for the given transfer function [32]. The problem is that the map from power spectrum to transfer functions is non-unique, and each factorization could lead to system realizations of differing dimension. For this reason, the concept of global minimality was introduced in [39] to select the transfer function with smallest system dimension. This raises the following question: Is global minimality sufficient to uniquely identify the transfer function from the power spectrum ? The answer is in general negative [56] , as discussed in [38,41] (see also Lemma 2 and Corollary 1 in [45] for a nice review). Our aim is to address these questions in the quantum case. In the following section we define an analogous notion of global minimality, and characterize globally minimal systems in terms of their stationary state. Afterwards we show that for SISO systems which admit a cascade realization the power spectrum and transfer function identification problems are equivalent.

A. Global minimality
As discussed earlier, in the time-dependent setting it is meaningful to restrict the attention to minimal systems, as they provide the lowest dimensional realizations which are consistent with a given input-output behavior. In the stationary setting however, it may happen that a minimal system can have the same power spectrum as a lower dimensional system. For instance if the input is the vacuum, and the system is passive then the stationary output is also vacuum and the power spectrum is trivial, i.e., the same as that of a zero-dimensional system. We therefore need to introduce a more restrictive minimality concept, as the stationary regime (power spectrum) counterpart of time-dependent (transfer function) minimality. The results of this section are valid for general MIMO systems and do not assume the existence of a cascade realization.
Definition 3. A system G = (S, C, Ω) is said to be globally minimal for input covariance V if there exists no lower dimensional system with the same power spectrum Ψ V . We call (G, V ) a globally minimal pair.
Before stating the main result of this section we briefly review some symplectic diagonalization results which will be used in the proof. Consider a k-modes cv system with canonical coordinatesc and a zero-mean Gaussian state with covariance matrix V := cc † . Any change of canonical coordinates which preserves the commutation relations is of the formc →c = Sc where S is a symplectic transformation S, cf. Definition 1. In the basisc , the state has covariance matrix V = SV S † . In particular there exists a symplectic transformation such that the modes c are independent of each other, and each of them is in a vacuum or a thermal state i.e. V i := c ic † i = ni+1 0 0 ni where n i is the mean photon number. We callc a canonical basis, and the elements of the ordered sequence n 1 ≤ · · · ≤ n k the symplectic eigenvalues of V . The latter give information about the state's purity: if all n i = 0 the state is pure, if all n i > 0 the state is fully mixed. More generally, we can separate the pure and mixed modes and write c = (c T p , c T m ) T . This procedure can be applied to the m input modes b, with covariance V (N, M ). Since the input is assumed to be pure, we have S in V (N, M )S † in = V vac where S in is a symplectic transformation and V vac is the vacuum covariance matrix. The interpretation is that any pure squeezed state looks like the vacuum when an appropriate symplectic "change of basis" is performed on the original modes.
Similarly, we can apply the above procedure to the stationary state of the system. Its covariance matrix P is the solution of the Lyapunov equation By an appropriate symplectic transformation we can change to a canonical basisȃ = S sysȃ such that a T = (a T p , a T m ). The system matrices are now A = S sys AS sys , C = CS sys . Note that this transformation is of the form prescribed by Theorem 1, but the interpretation here is that we are dealing with the same system seen in a different basis, rather than a different system with the same transfer function.
By combining the two symplectic transformations we see that any linear system with pure input can be alternatively described as a system with vacuum input and a canonical basis of creation and annihilation operators.
The following theorem links global minimality with the purity of the stationary state of the system. (1) The system is globally minimal if and only if the (Gaussian) stationary state with covariance P satisfying the Lyapunov equation (24) is fully mixed.
(2) A non-globally minimal system is the series product of its restriction to the pure component and the mixed component.
(3) The reduction to the mixed component is globally minimal and has the same power spectrum as the original system.
Proof. Let us prove the result first in the case S = 1.
First, perform a change of system and field coordinates as described above, so that the input is in the vacuum state, while the system modes decompose into its "pure" and "mixed" parts a T = (a T p , a T m ). Note that this transformation will alter the coupling and Hamiltonian matrices accordingly, but we still denote them Ω and C to simplify notations. Therefore, in this basis the stationary state of the system is given by the covariance P = R+1 0 0 R , R = 0 0 0 Rm and satisfies the Lyapunov equation (24). ( =⇒ ) We show that if the system has a pure component, then it is globally reducible. Let us write A ± and C ± as block matrices according to the pure-mixed splitting so that the Lyapunov equation (24) can be seen as a system of 16 block matrix equations. Taking the (1,1) and (1,3) blocks, which correspond to the a p a † p and a p a p components of the stationary state, one obtains Since Therefore, using this fact in Eq. (26) gives A pp + = 0, hence Ω pp + = 0. These two tell us that the pure part contains only passive terms.
Consider now the (1, 2) and (2, 3) blocks, which correspond to the a p a † m and a m a p components of the stationary state. From this, we get Similarly, Eq. (28) implies that A pm + = 0. Let G p := (1, Ω pp , C p ) be the system consisting of the pure modes, with Ω pp = ∆(Ω pp − , 0) and C p = ∆(C p − , 0). Let G m := (1, Ω mm , C m ) be the system consisting of the mixed modes with Ω mm = ∆(Ω mm − , Ω mm + ) and C m = ∆(C m − , C m + ). We can now show that the original system is the series product (concatenation) of the pure and mixed restrictions Indeed, using the fact that C p + = Ω pp + = A pm − = A pm + = 0, one can check that the series product has required matrices [10] C series =C p +C m = C and Ω series =Ω pp +Ω mm + Im (C mCp ) where the "tilde" notation stands for block matrices where only one block is nonzero, e.g.,C p = (C p , 0), and Im X := (X − X )/2i. Now, let Ξ p,m (s) denote the transfer functions of G p,m ; since the transfer function of a series product is the product of the transfer functions, we have Ξ(s) = Ξ m (s) · Ξ p (s). Furthermore, since G p is passive and the input is vacuum, we have Ψ p which means that the original system was globally reducible (not minimal).
( ⇐= ) We now show that if the system's stationary state is fully mixed, then it is globally minimal. The key idea is that a sufficiently long block of output has a finite symplectic rank (number of modes in a mixed state in the canonical decomposition) equal to twice the dimension of the system. Therefore the dimension of a globally minimal system is "encoded" in the output. This is the linear dynamics analog of the fact that stationary outputs of finite dimensional systems (or translation invariant finitely correlated states) have rank equal to the square of the system dimension (or bond dimension) [36]. To understand this property consider the system (S) together with the output at a long time 2T , and split the output into two blocks: A corresponding to an initial time interval [0, T ] and B corresponding to [T, 2T ]. If the system starts in a pure Gaussian state, then the S + A + B state is also pure. By ergodicity, at time T the system's state is close to the stationary state with symplectic rank d m . At this point the system and output block A are in a pure state so by appealing to the "Gaussian Schmidt decomposition" [52] we find that the state of the block A has the same symplectic eigenvalues (and rank d m ) as that of the system. In the interval [T, 2T ] the output A is only shifted without changing its state, but the correlations between A and S decay. Therefore the joint S + A state is close to a product state and has symplectic rank 2d m . On the other hand we can apply the Schmidt decomposition argument to the pure bipartite system consisting of S + A and B to find that the symplectic rank of B is 2d m . By ergodicity, B is close to the stationary state in the limit of large times, which proves the assertion.
To extend the result to S = 1, instead perform the change of field co-ordinates V → S in S b V S in S b † at the beginning. The proof for this case then follows as above because in this basis S = 1.
This result enables one to check global minimality by computing the symplectic eigenvalues of the stationary state. If all eigenvalues are nonzero, then the state is fully mixed and the system is globally minimal. We emphasize that the argument relies crucially on the fact that the input is a pure state. For mixed input states and in particular classical inputs, the stationary state may be fully mixed while the system is non globally minimal. The next step is to find out which parameters of a globally minimal system can be identified from the power spectrum.

A. Power spectrum identifiability result
The main result of this section is the following theorem which shows that two globally minimal SISO systems have the same power-spectrum if and only if they have the same transfer function, and in particular are related by a symplectic transformation as described in Theorem 1.
Writing V as S 0 ( 1 0 0 0 ) S † 0 for some symplectic matrix S 0 , we express the power spectrum as S 0Ξi (s)V vacΞi (−s # ) † S † 0 , whereΞ(s) is the transfer function of the system 1, S 0 C, Ω and V vac is the vacuum input. As S 0 is assumed to be known, the original problem reduces to proving the same statement for systems with vacuum input. In this case the power spectrum is given by The transfer function is completely characterized by the elements in the top row of its matrix, i.e., Ξ − (s) and Ξ + (s). Also, Ξ − (s) and Ξ + (s) must be of the the form (21) and (22). Our first observation is that Ξ − (s) and Ξ + (s) in (21) and (22) cannot contain poles and zeros in the following arrangement: Ξ − (s) has a factor like and Ξ + (s) contains a factor like For if this were the case and assuming that this could be done k times, then our original system could be decomposed as a cascade (series product) of two systems.
(i) The first system is a k-mode passive system with transfer function where Ξ (1) .
Note that by Example 2 it is physical.
(ii) The second system has n − k modes and transfer function where Ξ (2) It can be shown that there exists a minimal physical quantum system with this transfer function (see Appendix B).
Since Ξ (1) (s) is passive, and hence this k-mode system is not visible from the power spectrum, while the power spectrum is the same as that of the lower dimensional system Ξ (2) (s). Therefore we have a contradiction to global minimality. We will now construct Ξ − (s) and Ξ + (s) directly from the power spectrum. This is equivalent to identifying their poles and zeros [57]. To do this we must identify all poles and zeros of Ξ − (s) and Ξ + (s) from the three quantities: First, all poles of Ξ − (s) and Ξ + (s) may be identified from the power spectrum. Indeed, due to stability, each pole in (34)-(36) can be assigned unambiguously to either Ξ − (s) or Ξ + (−s). However, cancellations between zeros and poles of the two terms in the product may lead to some transfer function poles not being identifiable, so we need to show that this is not possible. Suppose that a pole λ of Ξ − (s) is not visible from the power spectrum. This implies the following We consider two separate cases: λ nonreal or real. If λ is nonreal then from the symmetries of the poles and zeros in (21) and (22), Ξ − (s) will contain a term like and Ξ + (s) will contain a term like By the argument above, the system is not globally minimal as there will be a mode of the system that is not visible in the power spectrum. Therefore all nonreal poles of Ξ − (s) may be identified. A similar argument ensures that all poles of Ξ + (s) are visible in the power spectrum.
If λ is real, we show that Ξ − (s) and Ξ + (s) will have terms of the form (37) and (38) and the result will follow. Indeed since λ is a pole of Ξ − (s), the denominator of (21) must have a second root at λ since the first cancels with the term (s − λ) which comes together with (s + λ) in the numerator. But then, Ξ + (s) must also have a pole at λ since otherwise |Ξ − (−iω)| 2 − |Ξ + (−iω)| 2 = 1 could not hold. A similar argument holds for a real pole of Σ + .
Therefore we conclude that all poles of Ξ ± (s) can be identified from the power spectrum, and we focus next on the zeros. Unlike the case of poles, it is not clear whether a given zero in any of these plots belongs to the factor on the left or the factor on the right in each of these equation (34), etc]. Since the poles of Ξ − (s) and Ξ + (s) may be different due to cancellations in (21) and (22), it is convenient here to add in "fictitious" zeros into the plots (34)-(36) so that Ξ − (s) and Ξ + (s) have the same poles. Note that these fictitious poles and zeros would have been present in (21) and (22) before simplification. From this point onwards, the zeros in (34)-(36) will refer to this augmented list which includes the additional zeros.
Real zeros.
In general the real zeros of Ξ − (s) and Ξ + (s) come in pairs ±λ [see Eqs. (21), (22)], unless a pole and zero (or more than one) cancel on the negative real line. Our task here is to distinguish these two cases from plots (34)- (36). Ξ − (s) has either (i) zeros at ±λ, or (ii) a zero at λ > 0 but not at −λ.
In case (i) (34) will have a double zero at each ±λ, whereas in case (ii) (34) will have a single zero at ±λ. We need to be careful here in discriminating cases (i) and (ii) on the basis of the zeros of (34). For example, a double zero at λ in (34) could be a result of one case (i) or two case (ii) in Ξ − (s). More generally, we could have an nth order zero at λ and as a result even more degeneracy is possible. A similar problem arises for the zeros of Ξ + (s) in (36).
Our first observation here is that it is not possible for both Ξ − (s) and Ξ + (s) to have zeros at ±λ (taking λ > 0 without loss of generality). If this were possible then by using the symplectic condition |Ξ − (−iω)| 2 − |Ξ + (−iω)| 2 = 1 and the fact that we are assuming that Ξ − (s) and Ξ + (s) have the same poles tells us that Ξ − (s) and Ξ + (s) must both have had double poles at −λ. The upshot is that Ξ − (s) and Ξ + (s) will have terms of the form (30) and (31), which is a contradiction. Now, suppose (34) has n zeros at λ > 0 and (36) has m zeros at λ > 0. Then we know that Ξ − (s) must have n−p 2 zeros at −λ and n+p 2 zeros at λ. Also, Ξ + (s) must have m−q 2 zeros at −λ and m+q 2 zeros at λ. The goal here is to find p and q because if these are known then it is clear that there must be n−p 2 ( m−q 2 ) type (i) zeros and p (q) type (ii) zeros in Ξ − (s) (Ξ + (s)).
By the observation above it is clear that either p = n or q = m. Also, in (35) there will be n+m+p−q 2 zeros at λ and n+m+q−p 2 zeros at −λ. Hence q − p is known at this stage. Finally, it is fairly easy to convince ourselves that if p = n but one concludes that q = m (or vice versa) and using the value of q − p leads to a contradiction. Hence p and q can be determined uniquely. For example, if n = 2, m = 5, q = 2 and p = 3 so that q = n and q − p = −1. Then assuming wrongly that p = 5 and using q − p = −1 it follows that q = 4 and so n must be 6, which is incorrect.
Having successfully identified all real zeros, we now show how to identify the zeros of Ξ − (s) and Ξ + (s) away from the real axis.
Comparing the zeros of (34) with those of (35) we find two cases in which the zeros can be assigned directly.
(i) Case 1: Let z be a zero of (34) that is not a zero of (35). Then z must be a zero of Ξ − (−s # ) # . Hence −z # is a zero of Ξ − (s).
(ii) Case 2: Let w be a zero of (35) that is not a zero of (34). Then w must be a zero of Ξ + (−s) # . Hence −w is a zero of Ξ + (s).
The question now is whether this procedure enables one to identify all zeros. Suppose that there is a zero v that is common to both of these plots. Then −v # must also be a zero of (34). Now, if −v # is not a zero of (35) then v is identifiable as belonging to Ξ − (s).
Therefore we can restrict our attention to the case that the zero pair {v, −v # } is common to both plots. Note that in this instance the list of zeros of (36) will also contain {v, −v # }. Assume without loss of generality that v is in the right half complex plane. Note that there cannot be a second zero pair {u, −u # } such that u = v # . If this were the case then either {v, −v} will be zeros of Ξ − (s) and {−v # , v # } will be zeros of Ξ + (s), or {u, −u} will be zeros of Ξ − (s) and {−u # , u # } will be zeros of Ξ + (s). In either case by using the condition |Ξ − (−iω)| 2 − |Ξ + (−iω)| 2 = 1 for all ω and the fact that Ξ − (s) and Ξ + (s) have the same poles by assumption, it follows that Ξ − (s) and Ξ + (s) will have terms of the form (30) and (31), which contradicts global minimality. Finally, under the assumptions that the zero pair {v, −v # } is common to both (35) and (34) with no second pair at {u, −u # } such that u = v # , then we can conclude that v must be a zero of Ξ − (s). For if this were not the case and so −v # were a zero of Ξ − (s) then there must be another zero of Ξ − (s) at v # (since pole-zero cancellation cannot occur in the right-half plane). Also from (35) this would require that Ξ + (s) has a zero at −v (hence also v). Therefore we have a contradiction to the fact that there is no second pair at {u, −u # } such that u = v # .
Therefore we have successfully identified all zeros of the transfer function away from the real axis, which completes the proof.
The theorem says that if a SISO system is globally minimal then the power spectrum is as informative as the transfer function. The result also gives a constructive method to check global minimality. Further it enables one to construct the transfer function of the system's globally minimal part. From this, one can then construct a system realization of this globally minimal restriction, using the results from Sec. III B. We call this realization method indirect because one first finds a transfer function fitting the power spectrum before constructing the system realization. Note that the work here also extends a result in [23]. There, conditions were derived to determine when the stationary state of the linear system is pure. Here, by means of the previous theorem, we have established a test to determine if there is a subsystem with a pure stationary state.
Remark 3. For general input V = S 0 V vac S † 0 , clearly systems of the form (S 0 ∆ (C − , 0) , ∆ (Ω − , 0)) have trivial power spectrum. Theorem 3 says that these are the only such systems (up to symplectic equivalence in Theorem 1).

Remark 4.
We have assumed that the scattering or squeezing matrix, S, for a system is the identity in this result. In fact the scattering or squeezing matrix is not always identifiable from the power spectrum. For example, a zero mode system with a single scattering term S = e iπ 0 0 e −iπ will have trivial power spectrum.

B. Power spectrum identification of passive QLSs
In this section we consider the special case of a minimal passive SISO QLSs. As noted before, we can therefore drop the doubled-up notation, cf. Eqs (13) and (14). For simplicity we will denoted C := C − , Ω := Ω − , and choose S = 1 so that the transfer function is where A = −iΩ − 1 2 C † C and its spectrum is σ(A) := {λ 1 , . . . , λ n }. The transfer function is a monic rational function in s, with poles p i = λ i in the left half plane, and zeros z i = −p i # = −λ # i in the right half plane. If the input state is vacuum then the power spectrum is trivial (Φ V = V ) and the only globally minimal systems are the trivial ones (zero internal modes). For this reason we restrict our attention to squeezed inputs, i.e., M = 0 in the input covariance. (1) The following are equivalent: (i) the system is globally minimal; (ii) the stationary state of the system is fully mixed; (iii) A and A † have different spectra, i.e., σ(A) ∩ σ(A † ) = ∅; (iv) A does not have real, or pairs of complex conjugate eigenvalues.
(2) Let P be the set of all eigenvalues of A that are either real or come in complex-conjugate pairs. A globally minimal realization of the system is given by the series product of one mode systems G m,i = (c i = 2|Reλ i |, Ω i = −Imλ i ) for indices i such that λ i / ∈ P. Proof.
(1) For passive SISO systems the only nontrivial contribution to the power spectrum is from off-diagonal element, In the above expression, zero-pole cancellations occur if and only if σ(A) ∩ σ(A † ) = ∅, or equivalently if A has a real eigenvalue or a pair of complex conjugate eigenvalues.
If no zero-pole cancellations occur, then σ(A) can be identified from Ξ(s)Ξ(−s) and the transfer function can be reconstructed. In this case the system is globally minimal.
If cancellations do occur then this happens in one of the two types of situations: (a) real eigenvalue: if λ i ∈ R then the corresponding term in the above product cancels (b) complex conjugate pairs: if λ i = λ # j then the i and j terms in the product cancel against each other.
In both cases, the remaining power spectrum has the same form, and can be seen as the power spectrum of a series product of one-dimensional passive systems, with dimension smaller than n, and therefore the system is not minimal.
This shows the equivalence of (i), (iii) and (iv) while the equivalence of (i) and (ii) follows from Theorem 2.
(2) The discussion so far shows that the transfer function factorizes as the product Ξ(s) = Ξ m (s)Ξ p (s) of a part corresponding to eigenvalues λ i ∈ P, which has a trivial power spectrum due to zero-pole cancellations, and the part corresponding to the complement which does not exhibit any cancellations. A system with transfer function Ξ(s) can be realized as series product G m G p of two separate passive systems with transfer functions Ξ m (s) and Ξ p (s). As argued before, G p has a pure stationary state which is uncorrelated to G m or the output, while G m has a fully mixed state which is correlated to the output.
Since G p does not contribute to the power spectrum, a globally minimal realization is provided by G m , Each fraction in (39) represents a bona fide PQLS G m,i with Hamiltonian and coupling parameters Ω i = −Imλ i and 1/2|c i | 2 = −Reλ i .
For PQLSs we now see that it is possible to construct a globally minimal realization of the PQLS directly from the power spectrum. Moreover, global minimality of PQLSs may be completely understood in terms of the spectrum of the system matrix A, just as was the case for minimality, stability, observability, and controllability [1,21]. An immediate corollary of this is the following. (1) the input is vacuum; (2) the eigenvalues of A are real or come in complexconjugate pairs.
From Theorem 4 there are two types of "elementary" systems that are not identifiable from the power spectrum for arbitrary input V (N, M ). Written in the doubled-up notation, these are either: (i) one mode systems of the form G 1 = (∆(c, 0), 0), or (ii) two mode systems of the form G 2 = (∆(c, 0), ∆(Ω − , 0)) (∆(c, 0), ∆(−Ω − , 0)). Either way it is not immediately obvious whether these systems are consistent with the nonidentifiable systems in Theorem 3. As an example we will show that this is indeed the case in the case of G 1 (G 2 is similar).
Example 3. Consider system G 1 for input V (N, M ), which is known to have (trivial) power spectrum V (N, M ). Therefore, in the vacuum basis of the field the system will beG G 1 must be transfer function equivalent to the system (∆(c, 0), 0) in the vacuum basis. Therefore, because this system is passive we have consistency with Theorem (3).
Finally, it seems that the assumption of global minimality seems to be not very restrictive; we illustrate this in the form of an example.
Example 4. Consider the following SISO PQLS with two internal modes: , where x ∈ R. We examine for which values of x the system is globally minimal for squeezed inputs. One can first check that the system is minimal if and only if x = 4. In Fig. 2 we plot the imaginary parts of the eigenvalues of A and A † . By Theorem 4, the system is not globally minimal if any of the lines representing the eigenvalues of A intersect those of A † . There are four points of interest that have been highlighted in the figure: 1 x = 0: crossing of eigenvalues of A but not with eigenvalues of A † ; system is globally minimal.
In summary, there were only two values of x for which the system is not globally minimal.

C. Global minimality with entangled inputs
Here we show that using an additional ancillary channel with an appropriate design of input makes it possible to identify the transfer function from the power spectrum for all minimal systems.
Consider the setup in Fig. 3, where a pure entangled input state is fed into a SISO QLS and an additional ancillary channel. The 2 × 2 blocks of the input V (N, M ) are The doubled-up transfer function is given by Now calculating the (2, 1) and (1, 4) entries of the power spectrum using (23), we obtain: Equivalently we may write these in matrix form as Hence if we choose |N 2 | = |M 2 | we may identify the transfer function of our SISO system uniquely. For example, such a choice of input would be N = x1 and M = 0 y y 0 with x(x + 1) = |y| 2 (the purity assumption). As one can see there are no requirements on the actual QLS other than minimality. Note that in the case of passive systems we need only that N 2 or M 2 be different from zero.

Remark 5.
Recall from the previous subsections that the maximum amount of information we may obtain about a PQLS from the power spectrum without the use of ancilla is that of the restriction to its globally minimal subspace. However, we have seen here that it is possible to construct a globally minimal pair, and hence obtain the whole transfer function simply by embedding the system in a larger space. To be clear here, there is no contradiction because the transfer function we are attempting to identify is the one in Eq. (41) rather than the SISO system Ξ(s).

VI. CONCLUSION
We have considered the identifiability of linear system using two contrasting approaches: (1) Time-dependent input (or transfer function) identifiability and (2) stationary inputs (or power spectrum) identifiability. In the time-dependent approach we characterized the equivalence class of systems with the same input-output data in Theorem 1, thus generalizing the results of [1] to active systems. We then outlined a method to construct a (minimal and physical) realization of the system from the transfer function. In fact, all results here hold for MIMO systems. In the stationary input regime, Theorem 2 showed that global minimality is equivalent to the stationary state of the system being fully mixed. Moreover, for a fixed pure input generically the transfer function may be constructed uniquely from the power spectrum under global minimality. A method was also given for how to do this in Theorem 3. Restricting to passive systems we saw that global minimality can be completely understood simply by considering the system matrix, A.
In particular, the transfer function can be constructed uniquely from the power spectrum if and only if none of the eigenvalues of A are real nor come in complexconjugate pairs (assuming that the input is squeezed). Finally, by using an ancillary channel it was shown that it is possible to identify any QLS uniquely from the transfer function. There are several directions to extend this work. First, it is expected that all results found for the stationary input approach can also be extended to (i) MIMO systems and (ii) those systems beyond the generic ones considered within this paper. We intend to address this in a future publication. Given that we now understand what is identifiable, the next step is to understand how well parameters can be estimated. In the time-dependent approach this has been done for passive systems in [1,22] but no such work exists for active systems or in the stationary approach at all. Last, it would be interesting to consider these identifiability problems in the more realistic scenario of noisy QLSs. In a QLS noise may be modelled by the inclusion of additional channels that cannot be monitored. Understanding what can be identified here will likely be far more challenging.

Appendix A: Finding a minimal classical realization
In this appendix a set of (nonphysical) minimal and doubled-up matrices (A 0 , B 0 , C 0 ) are found that realizes the transfer function (16), which describes a (minimal) physical system (A, C).
We assume that the matrix A for the n-mode minimal system, (A, C), possesses 2n distinct eigenvalues each with a nonzero imaginary part. This requirement can be seen to be generic in the space of all quantum systems [18]. Moreover, it can also be shown that if λ i is a complex eigenvalue of A with right eigenvector Ri Si and left eigenvector (U i , V i ), then λ # i is also an eigenvalue with right eigenvector U i , V i ∈ C n×1 and Σ n := 0n 1n 1n 0n . That is, for each eigenvalue and eigenvector, there exists a corresponding mirror pair. This property follows from the fact that A has the doubled-up form A := ∆ (A − , A + ).
We now construct a minimal realization called Gilbert's realization [43]. The only thing that we need to take care of is that the realization we obtain is of the doubled-up form.
As the transfer function may be written as .
we can perform a partial fraction expansion, so that As we show below, the matrices P i , Q i are rank 1. Therefore there exist matrices B i ∈ C 1×2 , B i ∈ C 1×2 , C i ∈ C 2×1 , and C i ∈ C 2×1 such that The Gilbert realization A 0 , B 0 , C 0 is A 0 := diag λ 1 , . . . , λ n , λ # 1 , . . . , λ # n ,  From the expression of the physical transfer function we have where W i are the rank-one matrices Having fixed B i and C i the matrices B i and C i can then be chosen as and so the matrices (A 0 , B 0 , C 0 ) are of the doubled-up type.
Note that using Gilbert's realization on MIMO systems can also be seen to give a minimal doubled-up realization, but we do not discuss this any further here.