Effects of parametric uncertainties in cascaded open quantum harmonic oscillators and robust generation of Gaussian invariant states

TL;DR

This paper analyzes how parametric uncertainties affect the generation of Gaussian invariant states in cascaded open quantum harmonic oscillators and proposes a method to optimize robustness against such uncertainties.

Contribution

It introduces an infinitesimal perturbation analysis for covariance matrices and develops a quadratic optimization approach for balancing oscillator realizations to enhance robustness.

Findings

Sensitivity of the purity functional can be minimized through symplectic transformations.

Effective solution for cascaded one-mode oscillators demonstrated numerically.

Connection established between sensitivity index and classical statistical distances.

Abstract

This paper is concerned with the generation of Gaussian invariant states in cascades of open quantum harmonic oscillators governed by linear quantum stochastic differential equations. We carry out infinitesimal perturbation analysis of the covariance matrix for the invariant Gaussian state of such a system and the related purity functional subject to inaccuracies in the energy and coupling matrices of the subsystems. This leads to the problem of balancing the state-space realizations of the component oscillators through symplectic similarity transformations in order to minimize the mean square sensitivity of the purity functional to small random perturbations of the parameters. This results in a quadratic optimization problem with an effective solution in the case of cascaded one-mode oscillators, which is demonstrated by a numerical example. We also discuss a connection of the sensitivity index with classical statistical distances and outline infinitesimal perturbation analysis for translation invariant cascades of identical oscillators. The findings of the paper are applicable to robust state generation in quantum stochastic networks.

Tables

Table 1. Table 1: Covariance bounds ( 213 ) for parametric uncertainties and contributions ( 216 ) to the mean square sensitivity index from three oscillators.

$k$	1	2	3
$a_k$	0.0222	0.0283	0.0067
$b_k$	0.0351	0.2898	0.0388
${\Psi }_k\left(I_2\right)$	37.9918	35.0268	19.5730
${\Psi }_k\left(S_k\right)$	34.6230	12.8844	14.4265
${\Psi }_k\left(S_k\right)/{\Psi }_k\left(I_2\right)$	0.9113	0.3678	0.7371

Equations

$$X:=\left[\begin{array}{c}{x}_1\\ ⋮\\ x_n\end{array}\right]$$

$$\mathrm{d}X=\mathcal{G}(X)\mathrm{d}t-i[X,L^{\mathrm{T}}]\mathrm{d}W$$

$$W:=\left[\begin{array}{c}{w}_1\\ ⋮\\ w_m\end{array}\right]$$

$$\mathrm{d}W\mathrm{d}{W}^{\mathrm{T}}:=(\mathrm{d}{w}_j\mathrm{d}{w}_k{)}_{1⩽j,k⩽m}=\Omega \mathrm{d}t,$$

$$\Omega$$

$$J$$

$$\mathbf{J}$$

$$[W(s),W(t{)}^{\mathrm{T}}]$$

$$=2i\min (s,t)J,s,t⩾0,$$

$$L:=\left[\begin{array}{c}{\ell }_1\\ ⋮\\ {\ell }_m\end{array}\right]$$

$$[X,L^{\mathrm{T}}]:=([x_j,{\ell }_k]{)}_{1⩽j⩽n,1⩽k⩽m}$$

$$\mathcal{G}(\xi):=i[H,\xi]+\frac{1}{2}\sum_{j,k=1}^{m}\omega_{jk}\big{(}[\ell_{j},\xi]\ell_{k}+\ell_{j}[\xi,\ell_{k}]\big{)},$$

$$\mathrm{d}\xi =\mathcal{G}(\xi )\mathrm{d}t-i[\xi ,L^{\mathrm{T}}]\mathrm{d}W.$$

$$\xi (t)=U(t{)}^{†}(\xi (0)\otimes {\mathcal{I}}_{\mathfrak{F}})U(t)$$

$$\mathrm{d}U(t)=-U(t)\Big{(}i(H(t)\mathrm{d}t+L(t)^{\mathrm{T}}\mathrm{d}W(t))+\frac{1}{2}L(t)^{\mathrm{T}}\Omega L(t)\mathrm{d}t\Big{)},$$

$$X(t)=U(t{)}^{†}(X(0)\otimes {\mathcal{I}}_{\mathfrak{F}})U(t).$$

$$Y(t):=\left[\begin{array}{c}{y}_1(t)\\ ⋮\\ y_m(t)\end{array}\right]=U(t{)}^{†}({\mathcal{I}}_{{\mathfrak{H}}_0}\otimes W(t))U(t)$$

$$\mathrm{d}Y=2JL\mathrm{d}t+\mathrm{d}W$$

$$[X(t),Y(s{)}^{\mathrm{T}}]=0,t⩾s⩾0$$

$$[Y(s),Y(t{)}^{\mathrm{T}}]=2i\min (s,t)J,s,t⩾0.$$

$${\mathcal{W}}_u{\mathcal{W}}_v={\mathrm{e}}^{i{v}^{\mathrm{T}}\Theta u}{\mathcal{W}}_{u+v}$$

$${\mathcal{W}}_u:={\mathrm{e}}^{i{u}^{\mathrm{T}}X}={\mathcal{W}}_{-u}^{†},$$

$$[{\mathcal{W}}_u,{\mathcal{W}}_v]=-2i\sin (u^{\mathrm{T}}\Theta v){\mathcal{W}}_{u+v},$$

$$[X,X^{\mathrm{T}}]=2i\Theta .$$

$$\Theta :=\frac{1}{2}\mathbf{J}\otimes I_{n/2}=\frac{1}{2}\left[\begin{array}{cc}0& I_{n/2}\\ -I_{n/2}& 0\end{array}\right],$$

$$H$$

$$L$$

$$\mathrm{d}X$$

$$\mathrm{d}Y$$

$$A$$

$$B$$

$$C$$

$$A\Theta +\Theta A^{\mathrm{T}}+BJ{B}^{\mathrm{T}}$$

$$\Theta C^{\mathrm{T}}+BJ$$

Full text

1 Introduction

The present paper is concerned with robustness of state generation in a class of open quantum systems with respect to unavoidable uncertainties which accompany practical implementation of such systems. This issue is important in the context of the emerging quantum information and quantum computation technologies which exploit the potential resources of physical systems at atomic scales described by quantum mechanics. Compared to classical systems with real-valued state variables whose evolution obeys the laws of Newtonian mechanics, quantum systems have more complicated operator-valued dynamic variables which act on a Hilbert space and evolve according to unitary similarity transformations. This unitary evolution is specified by an operator-valued Hamiltonian which quantifies the self-energy of the system when it is isolated from the environment. Interaction with a classical measuring device modifies the internal state of the quantum system in a random fashion, which depends on the quantum observable being measured and makes noncommuting operators inaccessible to simultaneous measurement [15, 26, 47]. The inherently stochastic nature of quantum systems is reflected in their quantum probability theoretic description [3, 27, 37] which replaces scalar-valued classical probability measures with density operators (quantum states) acting on the same Hilbert space as the dynamic variables. Regardless of whether measurements are involved, this probabilistic description is more complicated in the case of interaction between several open quantum systems, especially if one of them is organised as an infinite reservoir of “elementary” systems representing a quantum field.

A unified theoretic framework for the modelling and analysis of open quantum systems interacting with external bosonic fields (such as nonclassical light) is provided by the Hudson-Parthasarathy quantum stochastic calculus [17, 37] (see also [14]). This approach represents the Heisenberg picture evolution of system operators in the form of quantum stochastic differential equations (QSDEs) driven by noncommutative counterparts of the classical Wiener process [20]. Reflecting quantized energy exchange, the quantum Wiener processes involve annihilation and creation operators acting on symmetric Fock spaces [39]. The energetics of the quantum system itself and its interaction with the external fields, which specifies the structure of the QSDEs, is captured by the system Hamiltonian, system-field coupling operators and the scattering matrix (which pertains to the photon exchange between the fields in multichannel settings). In combination with the theory of quantum feedback networks [10, 18], this approach allows for the modelling of a wide class of interconnected open quantum systems which interact with one another and the environment. In fact, coherent (measurement-free) quantum control and filtering for quantum systems by direct or field-mediated interconnection [19, 28, 29, 34, 43, 44, 61, 62, 66, 73] constitute a promising modern paradigm which can potentially outperform the traditional observation-actuation approach of classical control theory. The principal advantage of this paradigm is that it avoids the loss of quantum information which accompanies the conversion of operator-valued quantum variables into classical real-valued signals in the process of measurement [7, 69].

The interaction of an open quantum system with external fields can be arranged in a dissipative fashion so that the resulting quantum stochastic system has an invariant state with desired properties. This provides an alternative approach [70, 72] to the system state preparation which can otherwise be carried out by steering the system to the required state through varying the parameters of the Hamiltonian in an open-loop fashion or by using feedback [6], similarly to the classical terminal state control problem [45, 58]. Since quantum state generation via dissipation does not involve measurement-based feedback, it is more aligned with the coherent quantum control paradigm mentioned above. State generation is important for quantum computation protocols [30], which often need the quantum system of interest to be initialized in a certain class of states.

Engineering a quantum state with required properties is relevant both for finite-level systems (such as qubit registers) and quantum systems with continuous variables which find applications in quantum optical platforms of quantum computing and quantum information processing [30, 68]. An important class of such systems is constituted by open quantum harmonic oscillators (OQHOs) whose dynamic variables satisfy the canonical commutation relations (CCRs), similar to those of the quantum mechanical positions and momenta [26, 47], and are governed by linear QSDEs [42]. The linearity of these QSDEs comes from the CCRs and the fact that the Hamiltonian and coupling operators of the OQHO are quadratic and linear functions of the system variables. The dynamics of such oscillators interacting with bosonic fields in the vacuum state are, in many respects, similar to those of classical Gaussian Markov diffusion processes (including the Ornstein-Uhlenbeck process [20]) and are particularly tractable at the level of the first two moments of the system variables. For example, the linear dynamics preserve the Gaussian nature [38, 40] of the system state in time, provided the OQHO is initialized in such a state. Moreover, irrespective of whether the initial state is Gaussian, linear QSDE with a Hurwitz matrix ensures the weak convergence [2, 5] of the system state to a unique invariant Gaussian state. It is this property of OQHOs that enables them to be employed for generating Gaussian quantum states by allowing the system to evolve over a sufficiently long period of time [23, 70]. The dissipation, which is built in this state generation procedure, secures stability of the invariant state being achieved in the long run [36].

However, the practical realization of OQHOs, which involves quantum optical components such as cavities, beam splitters and phase shifters, is accompanied by modelling inaccuracies and implementation errors. Even if the energetics of such a system remains linear-quadratic, the energy and coupling matrices can deviate from their nominal values. The resulting Gaussian invariant states can have covariance matrices which differ from the theoretical predictions. These deviations can lead to a deterioration or loss of important properties of quantum states such as purity [55], which can be critical due to the role of the pure states as extreme points of the convex set of density operators. This gives rise to the issue of robustness of the Gaussian state generation in linear quantum stochastic systems.

This circle of problems is the main theme of the present paper, which is concerned with Gaussian invariant states for cascaded OQHOs driven by vacuum fields. Such chain-like networks of linear quantum stochastic systems are relatively easy to implement and find applications, for example, in the generation of pure states [23, 70]. Since the transfer functions of the component oscillators satisfy the physical realizability conditions [19, 49], their $\mathcal{H}_{\infty}$-norms are not less than one. This can potentially enhance the propagation of perturbations in the energy and coupling matrices of the subsystems over long cascades. To this end, we carry out an infinitesimal perturbation analysis of the covariance matrix for the invariant Gaussian state and the related purity functional subject to inaccuracies in the energy and coupling matrices of the component oscillators. Since these quantities depend in a “spatially causal” fashion on the energy and coupling matrices of the oscillators along the cascade, their Frechet derivatives with respect to the matrix-valued parameters are amenable to recursive computation. This computation employs a connection with conditional covariance matrices of auxiliary classical Gaussian random vectors (which are updated similarly to the covariance equations of the discrete-time Kalman filter) and variational techniques for solutions of algebraic Lyapunov and Sylvester equations [56, 61].

By modelling the parametric uncertainties in the component oscillators as small zero-mean classical random elements which are uncorrelated for different subsystems, the corresponding perturbations in the purity functional are quantified by a mean square sensitivity index which depends on a particular state-space realization of the system. This leads to the problem of balancing the subsystems (for example, one-mode oscillators) by using CCR-preserving symplectic similarity transformations of their realizations in order to minimize the sensitivity of the purity functional to such uncertainties. By using group theoretic properties of Lyapunov equations, we obtain a closed-form dependence of the mean square sensitivity index on the transformation matrices, which allows its minimization to be decomposed into independent optimization problems for the subsystems. Each of these lower dimensional problems is organised as the minimization of a quartic polynomial of the corresponding symplectic matrix, which, in the one-mode case, reduces to a quadratic optimization problem with an effective solution. We demonstrate this result by a numerical example of balancing a cascade of one-mode oscillators.

The specific choice of the optimality criterion for balancing the cascaded oscillators is not unique and can be based on different cost functionals. To this end, we also discuss a connection of the above criterion with the classical Fisher information distance [48] applied to Gaussian quantum states. Note that balanced realizations, developed previously for classical linear systems and their quantum counterparts [33], were mainly concerned with equating the controllability and observability Gramians, which comes from the Kalman duality principle. However, the quantum setting of the present paper employs different criteria pertaining to the infinitesimal perturbation analysis of the state generation. In the case of translation invariant cascades of identical oscillators, the propagation of perturbations in subsystems over such cascades is particularly amenable to analysis using the technique of spatial $z$-transforms [63] (see also [53, 54]) which is also considered in the present study. The results of this paper may also find applications to perturbation analysis and robust Gaussian state generation in linear quantum stochastic networks with more complicated architectures.

The paper is organised as follows. Section 2 provides a background material on quantum stochastic systems. Section 3 describes the class of cascaded linear quantum stochastic systems under consideration. Section 4 specifies the invariant Gaussian quantum state for the composite system and the purity functional. Section 5 provides recurrence equations for the steady-state covariances using the cascade structure of the system. Section 6 carries out an infinitesimal perturbation analysis of the purity functional with respect to the energy and coupling matrices. Section 7 describes a recursive computation of the appropriate Frechet derivatives using the cascade structure of the system. Section 8 quantifies mean square sensitivity of the purity functional with respect to random implementation errors as an optimality criterion and splits its minimization into independent problems. Section 9 provides an illustrative numerical example of balancing a cascade of one-mode oscillators. Section 10 makes concluding remarks. The appendices provide an additional material for completeness. Appendix A discusses a connection between the sensitivity index for the purity functional and the Fisher information distance applied to Gaussian states. Appendix B carries out an infinitesimal perturbation analysis of steady-state covariances for translation invariant cascades of identical oscillators.

2 Linear quantum stochastic systems

For the purposes of the subsequent sections, we will provide a background on a class of quantum stochastic systems, including open quantum harmonic oscillators which play the role of building blocks in linear quantum systems theory [42]. The evolution of such a system in time $t\geqslant 0$ is described in terms of an even number $n$ of dynamic variables $x_{1}(t),\ldots,x_{n}(t)$ assembled into a vector

[TABLE]

(vectors are organised as columns unless indicated otherwise, and the time arguments are often omitted for brevity). These system variables are time-varying self-adjoint operators111there also are alternative formulations which employ non-Hermitian operators such as the annihilation and creation operators. on a complex separable Hilbert space $\mathfrak{H}$ (whose structure is clarified below), with $x_{1}(0),\ldots,x_{n}(0)$ acting on a subspace $\mathfrak{H}_{0}$ which is referred to as the initial system space. Their dynamics are governed by a Markovian Hudson-Parthasarathy quantum stochastic differential equation (QSDE) [17, 37]

[TABLE]

whose structure is described below. Although it formally resembles classical SDEs [20], the QSDE (2) is driven by a vector

[TABLE]

of an even number $m$ of self-adjoint operator-valued quantum Wiener processes $w_{1},\ldots,w_{m}$ acting on a symmetric Fock space $\mathfrak{F}$. These processes model the external bosonic fields [14, 17, 37], interacting with the system (see Fig. 1).

Accordingly, the above space $\mathfrak{H}$ is organised as the tensor-product Hilbert space $\mathfrak{H}:=\mathfrak{H}_{0}\otimes\mathfrak{F}$ which provides a common domain for the system and field operators. The quantum Wiener processes in (3) satisfy the quantum Ito relations

[TABLE]

where

[TABLE]

Here, $i:=\sqrt{-1}$ is the imaginary unit, $\otimes$ denotes the Kronecker product of matrices, $I_{r}$ is the identity matrix of order $r$ (which will often be omitted when it is clear from the context), and $\mathbf{J}$ spans the space of antisymmetric matrices of order $2$. Also, the transpose $(\cdot)^{\mathrm{T}}$ acts on matrices of operators as if their entries were scalars. In contrast to the identity diffusion matrix of the standard Wiener process, $\Omega$ in (5) is a complex positive semi-definite Hermitian matrix with an orthogonal antisymmetric imaginary part $\mathrm{Im}\Omega=J$ (so that $J^{2}=-I_{m}$). In view of (4), the quantum Wiener processes $w_{1},\ldots,w_{m}$ do not commute with each other. Furthermore, their two-point commutator matrix is given by

[TABLE]

where $[\alpha,\beta]:=\alpha\beta-\beta\alpha$ is the commutator of linear operators $\alpha$ and $\beta$. In accordance with $W$ representing the $m$-channel input field, the vector

[TABLE]

in (2) consists of system-field coupling operators $\ell_{1},\ldots,\ell_{m}$ which are also self-adjoint operators on the space $\mathfrak{H}$. Since the entries of the commutator matrix

[TABLE]

are skew-Hermitian operators, the dispersion $(n\times m)$-matrix $-i[X,L^{\mathrm{T}}]$ in (2) consists of self-adjoint operators on $\mathfrak{H}$. The $n$-dimensional drift vector $\mathcal{G}(X)$ of the QSDE (2) is obtained by the entrywise application of the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) generator [9, 21], which acts on a system operator $\xi$ (a function of the system variables) as

[TABLE]

where $\omega_{jk}$ are the entries of the quantum Ito matrix $\Omega$ from (5). Here, $H$ denotes the system Hamiltonian which is also a self-adjoint operator on $\mathfrak{H}$. The energy operators $H$ and $\ell_{1},\ldots,\ell_{m}$ are functions (for example, polynomials) of the system variables $x_{1},\ldots,x_{n}$ and inherit dependence on time from them. The GKSL superoperator $\mathcal{G}$ in (10) is a quantum analogue of the infinitesimal generators of classical Markov diffusion processes [20, 57] and specifies the drift of the QSDE

[TABLE]

The specific structure of the drift and diffusion terms in (11) (and its particular case (2)) comes from the evolution

[TABLE]

of the system operator $\xi$, with $\xi(0)$ acting on the initial system space $\mathfrak{H}_{0}$. Here, $(\cdot)^{\dagger}$ denotes the operator adjoint, and $U(t)$ is a unitary operator which acts on the system-field space $\mathfrak{H}$ and is governed by the QSDE

[TABLE]

where $U(0)=\mathcal{I}_{\mathfrak{H}}$ is the identity operator on $\mathfrak{H}$. The operator $U(t)$, which is associated with the system-field interaction over the time interval from [math] to $t$, is adapted in the sense that it acts effectively on the subspace $\mathfrak{H}_{0}\otimes\mathfrak{F}_{t}$, where $\{\mathfrak{F}_{t}:\,t\geqslant 0\}$ is the Fock space filtration.

The QSDE (11) can be obtained from (12) and (13) by using the quantum Ito formula [17, 37] in combination with (5), unitarity of $U$ and commutativity between the forward Ito increments $\mathrm{d}W(t)$ and adapted processes (including $U$) taken at time $s\leqslant t$. The corresponding quantum stochastic flow at time $t$ involves a unitary similarity transformation which acts on operators $\zeta$ on the space $\mathfrak{H}$ as $\zeta\mapsto U(t)^{\dagger}\zeta U(t)$ and applies entrywise to vectors of such operators. In accordance with (12), the system variables evolve as

[TABLE]

The system-field interaction, which drives the unitary operator $U(t)$ in (13), produces the output fields $y_{1},\ldots,y_{m}$ (see Fig. 1) which are also time-varying self-adjoint operators on the space $\mathfrak{H}$. The output field vector evolves as

[TABLE]

and satisfies the QSDE

[TABLE]

which is obtained similarly to (11). The QSDEs (2) and (16) describe a particular yet important scenario of quantum stochastic dynamics with the identity scattering matrix, which corresponds to the absence of photon exchange between the fields [17, 37]. Endowed with additional features (for example, more general scattering matrices and related gauge processes which affect the dynamics of the unitary operator $U$ in (13)), such QSDEs are employed in the formalism for modelling feedback networks of quantum systems which interact with each other and the external fields [10, 18]. Irrespective of a particular form of (13), the unitary similarity transformation in (14) and (15) preserves the commutativity between the system and output field variables in the sense that

[TABLE]

(that is, future system variables commute with the past output variables). At the same time, the output variables $y_{1},\ldots,y_{m}$ do not commute with each other and, in view of (8) and (16), inherit from the input fields the two-point commutator matrix:

[TABLE]

This noncommutativity makes the output fields inaccessible to simultaneous measurement [15, 26, 47]. However, they can be fed in a measurement-free fashion as an input to other open quantum systems. The resulting coherent field-mediated interconnection gives rise to fully quantum communication channels in the form of cascades of (possibly distant) quantum systems. In Section 3, we will consider such connections of linear quantum stochastic systems which are used, in particular, for the generation of certain classes of quantum states (see, for example, [23]).

The open quantum system, described above, is referred to as a linear quantum stochastic system, or an open quantum harmonic oscillator (OQHO) with $\frac{n}{2}$ modes, if its Hamiltonian and the coupling operators are quadratic and linear functions of the system variables, respectively, and the latter satisfy the following form of canonical commutation relations (CCRs). More precisely, the system variables $x_{1},\ldots,x_{n}$ of the OQHO satisfy the Weyl CCRs

[TABLE]

for all $u,v\in\mathbb{R}^{n}$. Here, $\Theta:=(\theta_{jk})_{1\leqslant j,k\leqslant n}$ is a real antisymmetric matrix of order $n$ (we denote the subspace of such matrices by $\mathbb{A}_{n}$) which is assumed to be nonsingular. The CCRs (18) are formulated in terms of the unitary Weyl operators

[TABLE]

where $u^{\mathrm{T}}X=\sum_{k=1}^{n}u_{k}x_{k}$ is a linear combination of the system variables assembled into the vector (1), with the coefficients comprising the vector $u:=(u_{k})_{1\leqslant k\leqslant n}\in\mathbb{R}^{n}$, so that $u^{\mathrm{T}}X$ is also a self-adjoint operator. The relations (18) imply that

[TABLE]

which leads to the Heisenberg infinitesimal form of the Weyl CCRs, specified on a dense subset of the space $\mathfrak{H}$ by the commutator matrix

[TABLE]

For example, in the case when the system variables are $q_{1},\ldots,q_{n/2},p_{1},\ldots,p_{n/2}$ and consist of conjugate quantum mechanical position $q_{k}$ and momentum $p_{k}=-i\partial_{q_{k}}$ operators (with an appropriately normalized Planck constant) [26, 47] which satisfy $[q_{j},p_{k}]=i\delta_{jk}$ for all $j,k=1,\ldots,\frac{n}{2}$, with $\delta_{jk}$ the Kronecker delta, the CCR matrix $\Theta$ takes the form

[TABLE]

where the matrix $\mathbf{J}$ is given by (7). Not restricted to this particular case, the system Hamiltonian $H$ of the OQHO is a quadratic function and the system-field coupling operators $\ell_{1},\ldots,\ell_{m}$ in (9) are linear functions of the system variables:

[TABLE]

Here, $R:=(r_{jk})_{1\leqslant j,k\leqslant n}$ is a real symmetric matrix of order $n$ (the subspace of such matrices is denoted by $\mathbb{S}_{n}$), and $M\in\mathbb{R}^{m\times n}$. The matrices $R$ and $M$ will be referred to as the energy and coupling matrices, respectively. Due to the CCRs (19) and the linear-quadratic energetics (21) and (22), the QSDEs (2) and (16) become linear with respect to the system variables:

[TABLE]

where the matrices $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, $C\in\mathbb{R}^{m\times n}$ are computed as

[TABLE]

In combination with the symmetry of $R$ and antisymmetry of $J$ and $\Theta$ in (6) and (19), the specific dependence of the matrices $A$, $B$, $C$ on the energy and coupling matrices $R$ and $M$ in (25)–(27) is equivalent to the physical realizability (PR) conditions [19, 49]:

[TABLE]

The equality (28) is related to the preservation of the CCRs (19) in time, while (29) corresponds to (17).

The linearity of the QSDEs (23) and (24) makes the OQHO a basic model in linear quantum control [19, 34, 42]. If the initial system variables have finite second moments (that is, $\mathbf{E}(X(0)^{\mathrm{T}}X(0))<+\infty$), the linear dynamics (2) preserve the mean square integrability in time and lead to finite limit values of the first and second moments

[TABLE]

provided the matrix $A$ in (25) is Hurwitz. Here, $\mathbf{E}\zeta:=\mathrm{Tr}(\rho\zeta)$ is the quantum expectation over the system-field density operator $\varrho:=\varpi\otimes\upsilon$, which is the tensor product of the initial system state $\varpi$ and the vacuum state $\upsilon$ of the input bosonic fields on the Fock space $\mathfrak{F}$ [37]. The matrix $P$ in (30) coincides with the infinite-horizon controllability Gramian

[TABLE]

of the matrix pair $(A,B)$, which is a unique solution of an algebraic Lyapunov equation (ALE) due to the matrix $A$ being Hurwitz:

[TABLE]

Despite this connection with classical linear systems theory, the quantum covariance matrix in (30) satisfies $P+i\Theta\succcurlyeq 0$, which reflects the generalized Heisenberg uncertainty principle [15] and is a stronger property than the positive semi-definiteness of $P$ alone in the classical case. This property can also be obtained directly by combining the ALE (32) with the PR condition (28), which leads to the ALE $A(P+i\Theta)+(P+i\Theta)A^{\mathrm{T}}+B\Omega B^{\mathrm{T}}=0$ whose solution satisfies

[TABLE]

in view of the positive semi-definiteness of the quantum Ito matrix $\Omega$ in (5). Application of a matrix-valued version of the Plancherel theorem leads to the frequency-domain representation of (33):

[TABLE]

where $(\cdot)^{*}:=(\overline{(\cdot)})^{\mathrm{T}}$ denotes the complex conjugate transpose. Here, $\mathbb{C}\ni s\mapsto F(s)\in\mathbb{C}^{n\times m}$ is a rational transfer function from $W$ to $X$ given by

[TABLE]

A related transfer function $G$ from $W$ to $Y$, which, together with $F$, is associated with the QSDEs (23) and (24) of the OQHO, is given by

[TABLE]

In view of (35), the poles of $F$ and $G$ belong to the spectrum of the matrix $A$ and hence, are entirely in the left half-plane $\mathrm{Re}s<0$ due to $A$ being Hurwitz. Similarly to classical linear stochastic systems, these transfer functions relate the Laplace transforms of the quantum processes $X$, $Y$, $W$ (considered in the right-half plane $\mathrm{Re}s>0$) as

[TABLE]

Due to the PR conditions (28) and (29), the transfer function $G$ in (36) is $(J,J)$-unitary [49] in the sense that

[TABLE]

for all $\lambda\in\mathbb{R}$. Since $G$ has an identity feedthrough matrix $\lim_{s\to\infty}G(s)=I_{m}$, its $\mathcal{H}_{\infty}$-norm (in the appropriate Hardy space) satisfies $\|G\|_{\infty}\geqslant 1$.

If the matrix $A$ is Hurwitz, then the reduced system state converges to an invariant Gaussian quantum state [38, 40] in the sense of appropriately modified weak convergence of probability measures [2]. This property holds irrespective of whether the initial system state is Gaussian or has finite second-order moments of the system variables (which is essential for (30)) and is equivalent to the point-wise convergence of the quasi-characteristic function (QCF) [5]:

[TABLE]

Here, the matrix $P$ is given by (31), and $\|u\|_{P}:=\sqrt{u^{\mathrm{T}}Pu}=|\sqrt{P}u|$ denotes the corresponding weighted Euclidean norm. The QCF on the right-hand side of (41) is identical to the characteristic function of the classical Gaussian distribution in $\mathbb{R}^{n}$ with zero mean vector and covariance matrix $P$. However, as mentioned above, the quantum nature of the setting manifests itself in the stronger property $P+i\Theta\succcurlyeq 0$ of the matrix $P$ in (33).

The convergence (41) can be used in order to generate a zero-mean Gaussian state with a given quantum covariance matrix $P+i\Theta$ as an invariant state of the OQHO described by (23) and (24); see, for example, [70]. To this end, the energy and coupling matrices $R$ and $M$ have to be chosen so that the matrix $A$ in (25) is Hurwitz and (31) is satisfied for a given admissible matrix $P$. Any particular choice of $R$ and $M$ does not affect the CCR matrix $\Theta$ in (19) or the Gaussian nature of the invariant state and only influences the matrix $P$. Indeed, the energetics of the OQHO remains linear-quadratic and the QSDE, driven by the vacuum input fields, remains linear.222A wider class of perturbations of the Hamiltonian and coupling operators, leading to nonlinear QSDEs and non-Gaussian invariant states, is considered, for example, in [52, 65, 67]. However, the matrix $P$ determines other important properties of Gaussian states such as purity [55]. In application to the invariant Gaussian system state, with a reduced density operator $r$ (which is also a positive semi-definite self-adjoint operator of unit trace), the purity is quantified by

[TABLE]

The Gaussian state is pure if and only if the inequality in (42) is an equality, in which case, $\det P$ achieves its minimum value $\det\Theta$ in view of the matrix inequality in (33).

In the case when the dimension $n$ of the OQHO is high, of practical interest is scalability in the generation of pure (or nearly pure) states. This can be achieved, for example, by using the cascade architecture [7, 23] which allows such an OQHO to be assembled from relatively simple components (such as one-mode oscillators). At the same time, the resulting large number of subsystems makes it important to secure robustness of the purity functional with respect to the cumulative effect of perturbations in the individual components.

3 Cascaded open quantum harmonic oscillators

Consider the field-mediated cascade connection of $N$ OQHOs in Fig. 2, which are driven by

the $m$-channel quantum Wiener process $W$ in (3) as described in Section 2. The oscillators have output fields $Y_{1},\ldots,Y_{N}$ of the same dimension [10, 23]. For every $k=1,\ldots,N$, the $k$th oscillator $O_{k}$ of the cascade is endowed with the initial space $\mathfrak{H}_{k}$ and a vector $X_{k}$ of $n_{k}$ system variables. In accordance with (19), these system variables satisfy CCRs and commute with one another for different subsystems:

[TABLE]

for any time $t\geqslant 0$ and all $j,k=1,\ldots,N$, where $\Theta_{k}\in\mathbb{A}_{n_{k}}$ are nonsingular CCR matrices. The system and output field variables are governed by a set of $N$ coupled QSDEs

[TABLE]

where $Y_{0}$ is identified with $W$, as mentioned before. The matrices $A_{k}\in\mathbb{R}^{n_{k}\times n_{k}}$, $B_{k}\in\mathbb{R}^{n_{k}\times m}$, $C_{k}\in\mathbb{R}^{m\times n_{k}}$ of these QSDEs are given by

[TABLE]

similarly to (25)–(27). Here, $J\in\mathbb{A}_{m}$ is the imaginary part of the quantum Ito matrix $\Omega$ in (5) which is inherited by the output fields $Y_{1},\ldots,Y_{N}$ in (45) from the input quantum Wiener process $W$ in (3):

[TABLE]

for all $j,k=1,\ldots,N$. Also, similarly to (21) and (22), the energy and coupling matrices $R_{k}\in\mathbb{S}_{n_{k}}$ and $M_{k}\in\mathbb{R}^{m\times n_{k}}$ in (46)–(48) specify the system Hamiltonian $H_{k}$ and the vector $L_{k}$ of $m$ operators of system-field coupling of the $k$th oscillator $O_{k}$ to the output $Y_{k-1}$ of the preceding oscillator $O_{k-1}$ in the cascade (see Fig. 2):

[TABLE]

In accordance with (28) and (29), the state-space matrices $A_{k}$, $B_{k}$, $C_{k}$ in (46)–(48) satisfy the PR conditions

[TABLE]

which pertain to the preservation of the CCRs (43) and the commutativity between the future system variables and the past output variables of the cascaded oscillators:

[TABLE]

The above cascade can be regarded as a composite OQHO with the augmented vector of $n:=n_{1}+\ldots+n_{N}$ system variables

[TABLE]

which satisfies a linear QSDE obtained by combining the QSDEs (44) and (45):

[TABLE]

The output of this composite system is the output $Y_{N}$ of the $N$th oscillator $O_{N}$ in the cascade. Also, the matrices $\mathcal{A}\in\mathbb{R}^{n\times n}$, $\mathcal{B}\in\mathbb{R}^{n\times m}$, $\mathcal{C}\in\mathbb{R}^{m\times n}$ in (55) and (56) are the last elements

[TABLE]

of a matrix sequence $(\mathcal{A}_{k},\mathcal{B}_{k},\mathcal{C}_{k})_{1\leqslant k\leqslant N}$ computed recursively as

[TABLE]

with the initial conditions being the matrices of the first oscillator $O_{1}$ in the cascade:

[TABLE]

see also [23]. In view of (46)–(48) and (58), the blocks $A_{jk}\in\mathbb{R}^{n_{j}\times n_{k}}$ of the matrix $\mathcal{A}$ are given by

[TABLE]

for all $j,k=1,\ldots,N$. The corresponding augmented version of the PR conditions (51) and (52) takes the form

[TABLE]

where

[TABLE]

is a block-diagonal CCR matrix of order $n$ for the augmented vector $\mathcal{X}$ in (54) which represents the CCRs (43) as

[TABLE]

Similarly to (51) and (52), the PR conditions (63) and (64) are equivalent to the preservation of the CCRs (66) and the commutativity (53), respectively.

In accordance with [31, Theorem 5.1] and [32, Lemma 3], a combination of (46)–(48) with (57)–(62) and (65) leads to the following energy and coupling matrices $\mathcal{R}\in\mathbb{S}_{n}$ and $\mathcal{M}\in\mathbb{R}^{m\times n}$ of the composite OQHO:

[TABLE]

where the blocks $R_{jk}\in\mathbb{R}^{n_{j}\times n_{k}}$ are computed for all $j,k=1,\ldots,N$ as

[TABLE]

The energy matrices $R_{1},\ldots,R_{N}$ of the component oscillators enter the energy matrix $\mathcal{R}$ of the composite system only through its diagonal part

[TABLE]

which takes values in a Hilbert space $\mathfrak{R}$, isomorphic to the Hilbert space $\mathbb{S}_{n_{1}}\times\ldots\times\mathbb{S}_{n_{N}}$ with the direct-sum inner product (associated with the Frobenius inner product of real matrices $\langle L,M\rangle:=\mathrm{Tr}(L^{\mathrm{T}}M)$ and the corresponding norm $\|M\|:=\sqrt{\langle M,M\rangle}$). The specific quadratic dependence of the energy matrix $\mathcal{R}$ on the coupling matrices $M_{1},\ldots,M_{N}$ in (67) and (68) is closely related to the block lower triangularity of the matrix $\mathcal{A}$ of the augmented QSDE (55) due to the cascade structure of the quantum system being considered. While the matrix $\mathcal{A}$, which can be represented as

[TABLE]

depends linearly on the energy matrices $R_{1},\ldots,R_{N}$ and quadratically on the coupling matrices $M_{1},\ldots,M_{N}$ of the constituent OQHOs, the dependence of the matrices $\mathcal{B}$ and $\mathcal{C}$ in (55) and (56) on the coupling matrices is linear, similarly to (47) and (48):

[TABLE]

The Hamiltonian $\frac{1}{2}\mathcal{X}^{\mathrm{T}}\mathcal{R}\mathcal{X}$ and the vector $\mathcal{M}\mathcal{X}$ of operators of coupling with the external input $W$, which are computed for the composite OQHO in terms of the component parameters according to (67) and (68), can also be obtained by using the quantum feedback network formalism [10, 18] and the individual energy operators in (50).

4 Gaussian invariant quantum state

If the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz, then so is the block lower triangular matrix $\mathcal{A}$ in (57). In this case, in accordance with Section 2, the composite OQHO has a unique invariant state which is Gaussian with zero mean vector and the quantum covariance matrix $\mathcal{P}+i\mathbfit{\Theta}\succcurlyeq 0$, whose real part

[TABLE]

is the controllability Gramian of the pair $(\mathcal{A},\mathcal{B})$, similarly to (31) and (32), and is found uniquely by solving the ALE

[TABLE]

Similarly to (33) and (34), the state-space and frequency-domain representations of the invariant quantum covariance matrix are given by

[TABLE]

Here, $\mathcal{F}$ is the $\mathbb{C}^{n\times m}$-valued rational transfer function from the input quantum Wiener process $W$ to the system variables in $\mathcal{X}$ computed as

[TABLE]

where

[TABLE]

are the transfer functions from $Y_{k-1}$ to $X_{k}$ and $Y_{k}$, respectively, associated with the QSDEs (44) and (45) for the $k$th OQHO. Similarly to (37) and (38), these transfer functions relate the Laplace transforms of $X_{k}$ and $Y_{k}$ to that of $W$ in (39) in the right-half plane $\mathrm{Re}s>0$ by

[TABLE]

In view of the PR conditions (51) and (52), each of the functions $G_{1},\ldots,G_{N}$ in (77) is $(J,J)$-unitary in the sense of (40). The products of the transfer matrices in (75) are closely related to the concatenation of quantum systems [10, 18] in the case of linear QSDEs being considered [71]. Since the transfer functions $G_{k}$ in (77) have an identity feedthrough matrix $I_{m}$, so also do their products. Hence, their $\mathcal{H}_{\infty}$-norms satisfy

[TABLE]

which also follows from the $(J,J)$-unitarity333and the fact that such matrices form a group of $G_{k}$, whereby no contraction can be guaranteed in (75). This can potentially facilitate the propagation of modelling errors over the cascade (especially from the first oscillators towards the end in long cascades), which will be discussed for the translation invariant case in Appendix B.

Now, a given admissible real part $\mathcal{P}$ of the quantum covariance matrix $\mathcal{P}+i\mathbfit{\Theta}\succcurlyeq 0$ for the invariant Gaussian state of the cascaded OQHOs in (43)–(45) can be achieved by choosing the energy and coupling matrices $R_{1},\ldots,R_{N}$ and $M_{1},\ldots,M_{N}$ so as to make the matrices $A_{1},\ldots,A_{N}$ in (46) Hurwitz and to satisfy (72); see, for example, [23, 70]. Inaccuracies in these matrices leave the CCRs (66) intact and can only perturb the matrix $\mathcal{P}$ without destroying the Gaussian nature of the invariant state. In application to the invariant Gaussian state of the composite OQHO (with the reduced density operator $r$), the purity functional (42) takes the form

[TABLE]

where only the denominator depends on the matrices $R_{1},\ldots,R_{N}$ and $M_{1},\ldots,M_{N}$ of the oscillators. Note that the purity functional (81) is invariant with respect to the transformations

[TABLE]

for arbitrary matrices $S_{k}\in\mathbb{R}^{n_{k}\times n_{k}}$, satisfying $S_{k}\Theta_{k}S_{k}^{\mathrm{T}}=\Theta_{k}$ and forming the symplectic group $\mathrm{Sp}(\Theta_{k})$ (which agrees with the usual definition of the symplectic group [11, 35, 60] up to the matrix transpose). With the transfer functions $G_{k}$ in (77) remaining unchanged, (82) corresponds to the symplectic transformations

[TABLE]

of the system variables of the component oscillators which preserve the CCRs (66) since the matrix

[TABLE]

belongs to $\mathrm{Sp}(\mathbfit{\Theta})$, where $\mathbfit{\Theta}$ is given by (65). From the property $\det S=1$ for any such matrix $S$, it follows that the corresponding transformation

[TABLE]

leaves $\det(S\mathcal{P}S^{\mathrm{T}})=\det\mathcal{P}$ unchanged, and hence, the purity functional in (81) is indeed invariant under such transformations. Since (81) can be expressed in terms of the log determinant of the matrix $\mathcal{P}$ as

[TABLE]

the minimization of the functional $V$ and its sensitivity to the perturbations in the energy and coupling matrices is relevant to the purity of Gaussian quantum states generated through cascaded oscillators [23, 32, 70].

5 Recursive computation of the steady-state covariances

The cascade structure of the composite quantum system under consideration (see Fig. 2) leads to “spatial causality” in the dependence of the matrix $\mathcal{P}$ in (72) on the energy and coupling matrices of the component oscillators. In accordance with the partitioning of the vector $\mathcal{X}$ in (54), the matrix $\mathcal{P}$ can be split into blocks $P_{jk}\in\mathbb{R}^{n_{j}\times n_{k}}$ as

[TABLE]

and regarded as the last element in the sequence of matrices $\mathcal{P}_{1},\ldots,\mathcal{P}_{N}$ given by

[TABLE]

where

[TABLE]

The matrix $\mathcal{P}_{k}$ in (88) is the real part of the quantum covariance matrix $\mathcal{P}_{k}+i\mathbfit{\Theta}_{k}$ of the reduced Gaussian invariant state for the oscillators $O_{1},\ldots,O_{k}$, where $\mathbfit{\Theta}_{k}:=\mathop{\mathrm{diag}}_{1\leqslant j\leqslant k}(\Theta_{j})$ is the CCR matrix of their system variables constituting the vectors $X_{1},\ldots,X_{k}$, which corresponds to (65). The above mentioned spatial causality is the property that $\mathcal{P}_{j}$ is independent of $R_{k}$ and $M_{k}$ for all $1\leqslant j<k\leqslant N$. This property is closely related to the following recursive computation of the matrix $\mathcal{P}$.

Lemma 1**.**

Suppose the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz. Then the matrix $\mathcal{P}$ in (87) can be found by recursively computing the matrices $\mathcal{P}_{1},\ldots,\mathcal{P}_{N}$ in (88) through solving the ALEs

[TABLE]

with respect to the matrices $Q_{k}$ and $P_{kk}$ in (89) for all $k=2,\ldots,N$. Here, the initial condition $\mathcal{P}_{1}=P_{11}$ is obtained by solving the ALE

[TABLE]

$\square$**

Proof.

For any $k=1,\ldots,N$, the matrix $\mathcal{P}_{k}$, which pertains to the steady-state covariance dynamics of the system variables for the first $k$ oscillators in the cascade, satisfies the ALE

[TABLE]

where the matrices $\mathcal{A}_{k}$ and $\mathcal{B}_{k}$ are given by (58) and (59). The ALE (92) is obtained by letting $k=1$ in (93). For any $k=2,\ldots,N$, the structure of the matrices $\mathcal{A}_{k}$ and $\mathcal{B}_{k}$ in (58), (59) and $\mathcal{P}_{k}$ in (88) leads to

[TABLE]

Therefore, since the left-hand side of (93) is a symmetric matrix, this ALE splits into three equations consisting of (90), (91) and

[TABLE]

The last equation reproduces (93) with the preceding value of $k$, thus completing the proof by a standard induction argument. $\blacksquare$

Note that for any given $k=1,\ldots,N$, the energy and coupling matrices $R_{k}$ and $M_{k}$ of the $k$th oscillator influence the matrix $\mathcal{P}$ through the $k$th block rows and the $k$th block columns of the matrices $\mathcal{A}$ and $\mathcal{B}\mathcal{B}^{\mathrm{T}}$ in the ALE (73). The corresponding “cross-shaped” fragments of $\mathcal{A}$ and $\mathcal{B}\mathcal{B}^{\mathrm{T}}$ are described by

[TABLE]

where use is made of (58) and (59) together with an auxiliary matrix

[TABLE]

At the level of covariances, there is a connection between the quantum system variables444which, as mentioned before, are not accessible to simultaneous measurement and conditional averaging because of their noncommutativity [15] and classical random variables. More precisely, the real parts of the steady-state covariances of the system variables in the cascaded OQHOs can be reproduced by considering a sequence of jointly Gaussian classical random vectors $\xi_{1},\ldots,\xi_{N}$ of dimensions $n_{1},\ldots,n_{N}$ in the form

[TABLE]

where $\omega$ is a standard Wiener process in $\mathbb{R}^{m}$, and the improper integral is convergent (in particular, in the mean square sense) since the matrix $\mathcal{A}_{k}$ is Hurwitz. Indeed, the covariance matrices of these auxiliary random vectors coincide with the blocks of the matrix $\mathcal{P}$ in (72) and (87)–(89):

[TABLE]

for all $j,k=1,\ldots,N$, in terms of which the functional $V$ in (86) admits the decomposition

[TABLE]

Here, the matrix $\Pi_{k}$ is the Schur complement [16] of the block $\mathcal{P}_{k-1}$ in the matrix $\mathcal{P}_{k}$ in (88), which is given by

[TABLE]

and coincides with the conditional covariance matrix of the random vector $\xi_{k}$ given the “past history” $\Xi_{k-1}$ (in the sense of the spatial parameter $k$). Also,

[TABLE]

is the corresponding predictor. In (98) and (99), use is also made of (95), (96) and the structure of conditional distributions for jointly Gaussian random vectors, which plays an important role in linear stochastic filtering [1, 22].

The above mentioned spatial causality manifests itself in the fact that, for all $1\leqslant j<k\leqslant N$, the Schur complement $\Pi_{j}$ in (98) is independent of the energy and coupling matrices $R_{k}$ and $M_{k}$ of the $k$th oscillator, and hence, so is the quantity $V_{j}$ in (97).

6 Infinitesimal perturbation analysis of the purity functional

As mentioned in Section 4, inaccuracies in the energy and coupling matrices of the cascaded OQHOs can affect the purity functional (81) for the Gaussian invariant state of the system. Its sensitivity to infinitely small perturbations can be described in terms of the Frechet derivatives of the functional $V$ in (86) with respect to the matrices $R_{1},\ldots,R_{N}$ and $M_{1},\ldots,M_{N}$ which constitute the matrices $\mathbfit{R}$ and $\mathcal{M}$ in (69) and (67). To this end, we will first establish a preliminary lemma which computes such derivatives for the invariant covariance matrix of the composite system. More precisely, the condition that the matrix $\mathcal{A}$ of the composite quantum system (55) is Hurwitz ensures smooth dependence of the matrix $\mathcal{P}$ in (72) on the $\mathfrak{R}\times\mathbb{R}^{m\times n}$-valued pair $(\mathbfit{R},\mathcal{M})$ comprising the energy and coupling matrices of the cascaded oscillators. This smoothness leads to a well-defined Frechet derivative

[TABLE]

which (for a given pair $(\mathbfit{R},\mathcal{M})$) is a linear operator acting from the space $\mathfrak{R}\times\mathbb{R}^{m\times n}$ to $\mathbb{S}_{n}$. The first variation of the matrix $\mathcal{P}$ is expressed in terms of $\Lambda$ as

[TABLE]

where $\partial_{\mathbfit{R}}\mathcal{P}$ and $\partial_{\mathcal{M}}\mathcal{P}$ denote the corresponding partial Frechet derivatives. For what follows, we associate with arbitrary Hurwitz matrices $\alpha$ and $\beta$ a linear operator $\mathbf{L}_{\alpha,\beta}$ which maps an appropriately dimensioned matrix $\gamma$ to

[TABLE]

which is a unique solution $\sigma$ of an algebraic Sylvester equation (ASE), or a generalized ALE $\alpha\sigma+\sigma\beta^{\mathrm{T}}+\gamma=0$; see, for example, [8]. In the case when $\alpha=\beta$, which corresponds to standard ALEs, an abbreviated notation

[TABLE]

will be used. Also, we denote by $[\![\![\alpha,\beta]\!]\!]$ a “sandwich” operator which acts on an appropriately dimensioned matrix $\gamma$ as

[TABLE]

Furthermore, it will be convenient to denote the symmetrizer and antisymmetrizer of square matrices by

[TABLE]

Also, we will need an auxiliary linear operator $\Gamma:=\Gamma_{\mathcal{M}}$, which is associated with the coupling matrix $\mathcal{M}$ in (67) and maps its variation $\delta\mathcal{M}$ to a block lower triangular matrix $\Gamma(\delta\mathcal{M}):=(\Gamma_{jk}(\delta\mathcal{M}))_{1\leqslant j,k\leqslant N}$ whose blocks are given by

[TABLE]

The first variation of the block $A_{jk}$ of the matrix $\mathcal{A}$ in (62) with respect to the matrix $\mathcal{M}$ is expressed in terms of (106) as

[TABLE]

for all $j,k=1,\ldots,N$. Therefore, the operator $\Gamma$ allows the partial Frechet derivative of $\mathcal{A}$ with respect to $\mathcal{M}$ to be represented as the composition

[TABLE]

where (65) and (104) are used. The following lemma carries out an infinitesimal perturbation analysis for the steady-state covariances of the quantum system.

Lemma 2**.**

Suppose the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz. Then the partial Frechet derivatives $\partial_{\mathbfit{R}}\mathcal{P}$ and $\partial_{\mathcal{M}}\mathcal{P}$ in (101) are computed as

[TABLE]

in terms of the operators (102)–(107). $\square$

Proof.

The proof is carried out by using the algebraic techniques of [56, 61] justified by the smooth dependence of $\mathcal{P}$ on the energy and coupling matrices under the condition that $\mathcal{A}$ is Hurwitz. More precisely, the first variation of the ALE (73) leads to

[TABLE]

Since the matrix $\mathcal{P}$ is symmetric, (110) can be represented as

[TABLE]

where use is made of the symmetrizer $\mathbf{S}$ from (105). The relation (111) is an ALE with respect to the matrix $\delta\mathcal{P}$, whose solution can be expressed as

[TABLE]

in terms of the operator $\mathbf{L}_{\mathcal{A}}$ given by (102) and (103). In view of (62) and (70), the first variation of the matrix $\mathcal{A}$, as a function of the matrices $\mathbfit{R}$ and $\mathcal{M}$ from (69) and (67), is computed as

[TABLE]

where use is made of the operator $\Gamma$ from (106) and (107). By a similar reasoning, the first variation of the matrix $\mathcal{B}$ in (71) with respect to $\mathcal{M}$ is given by

[TABLE]

Substitution of (113) and (114) into (112) relates $\delta\mathcal{P}$ to $\delta\mathbfit{R}$ and $\delta\mathcal{M}$ by

[TABLE]

where use is also made of the antisymmetry of the CCR matrix $\mathbfit{\Theta}$ in (65). By comparing (115) with (101), it follows that the partial Frechet derivatives $\mathcal{P}$ with respect to $\mathbfit{R}$ and $\mathcal{M}$ take the form (108) and (109). $\blacksquare$

The proof of Lemma 2 leads to the following group theoretic property of the linear operator $\Lambda$ in (100) with respect to the symplectic similarity transformations of the cascaded oscillators described by (82) and (83).

Lemma 3**.**

Suppose the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz. Then, for any symplectic matrix $S\in\mathrm{Sp}(\mathbfit{\Theta})$ in (84), the Frechet derivatives $\Lambda$ in (100), evaluated at the original and transformed pairs $(\mathbfit{R},\mathcal{M})$ and $(S^{-\mathrm{T}}\mathbfit{R}S^{-1},\mathcal{M}S^{-1})$, are related by

[TABLE]

$\square$**

Proof.

The property (116) is inherited by $\Lambda$ from solutions of the ALE (110). More precisely, if the matrices $\mathbfit{R}$ and $\mathcal{M}$ and their variations $\delta\mathbfit{R}$ and $\delta\mathcal{M}$ are transformed as

[TABLE]

then the corresponding matrices $\mathcal{A}$, $\mathcal{B}$, $\mathcal{P}$ and their first variations $\delta\mathcal{A}$, $\delta\mathcal{B}$, $\delta\mathcal{P}$ in (113), (114), (110) are transformed as

[TABLE]

due to the symplectic property of the matrix $S\in\mathrm{Sp}(\mathbfit{\Theta})$. A combination of (117)–(121) with (101) implies that the operator $\Lambda$ in (100) satisfies

[TABLE]

The relation (116) now follows from (122) by applying an appropriate inverse transformation to the variations $\delta\mathbfit{R}$ and $\delta\mathcal{M}$ of the independent variables. An alternative (and somewhat less intuitive) way to establish (116) is to use the operator identities

[TABLE]

which hold for any matrix $S\in\mathrm{Sp}(\mathbfit{\Theta})$ in (84). Indeed, by applying (123)–(127) to the right-hand sides of (108) and (109), it follows that the operators $\partial_{\mathbfit{R}}\mathcal{P}$ and $\partial_{\mathcal{M}}\mathcal{P}$ are transformed as

[TABLE]

which is an equivalent form of (116). $\blacksquare$

Lemma 3 simplifies the “recalculation” of the Frechet derivative $\Lambda$ under the symplectic similarity transformations of a particular realization of the cascaded oscillators. Indeed, the linear operator $\Lambda$ is modified by applying appropriate linear transformations $\mathfrak{S}_{S^{-1}}$ and $\mathfrak{P}_{S}$ to its domain $\mathfrak{R}\times\mathbb{R}^{m\times n}$ and range $\mathbb{S}_{n}$ as

[TABLE]

which is merely a concise representation of (128) and (129) or their equivalent form (116). Here, the transformations

[TABLE]

are parameterized by the symplectic matrix $S\in\mathrm{Sp}(\mathbfit{\Theta})$ in (84) and describe group homomorphisms (whereby $\mathfrak{S}_{S^{-1}}=\mathfrak{S}_{S}^{-1}$).

The following theorem, which is concerned with infinitesimal perturbation analysis of the purity functional (81), applies Lemma 3 in order to obtain the corresponding transformations for the partial Frechet derivatives of $V$ in (86) with respect to the energy and coupling matrices $R_{1},\ldots,R_{N}$ and $M_{1},\ldots,M_{N}$ of the component oscillators:

[TABLE]

Up to a constant factor of $-\frac{1}{2}$, the matrices $\rho_{k}$ and $\mu_{k}$ describe the corresponding logarithmic Frechet derivatives of the purity functional. Before actually computing $\rho_{k}$ and $\mu_{k}$ (which are fairly complicated rational functions of the entries of the energy and coupling matrices), the theorem shows that, under the symplectic similarity transformations, these derivatives are at most quadratic polynomials of the transformation matrices.

Theorem 4**.**

Suppose the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz. Then, for any symplectic matrix $S\in\mathrm{Sp}(\mathbfit{\Theta})$ in (84), the Frechet derivatives of the functional $V$ with respect to the energy and coupling matrices $R_{1},\ldots,R_{N}$ and $M_{1},\ldots,M_{N}$ in (130) are transformed as

[TABLE]

for all $k=1,\ldots,N$. $\square$

Proof.

We will combine Lemma 3 with the following Frechet derivative on the Hilbert space $\mathbb{S}_{n}$ (see, for example, [13, 16]):

[TABLE]

This allows the first variation of the functional $V$ in (86) (as a composite function $(\mathbfit{R},\mathcal{M})\mapsto\mathcal{P}\mapsto\ln\det\mathcal{P}$) to be computed as

[TABLE]

whence

[TABLE]

in view of (130). From the transformations (128) and (129), it follows that the corresponding adjoint operators $\partial_{\mathbfit{R}}\mathcal{P}^{\dagger}$ and $\partial_{\mathcal{M}}\mathcal{P}^{\dagger}$ are modified as

[TABLE]

where use is also made of the relation

[TABLE]

for the operator (104) which yields $[\![\![S,S^{\mathrm{T}}]\!]\!]^{\dagger}=[\![\![S^{\mathrm{T}},S]\!]\!]$ and $[\![\![I_{m},S]\!]\!]^{\dagger}=[\![\![I_{m},S^{\mathrm{T}}]\!]\!]$. A combination of (135) and (136) with the transformation (85) implies that the Frechet derivatives in (133) and (134) are transformed as

[TABLE]

The relations (131) can now be obtained as a block-wise form of (138) and (139) due to the block diagonal structure of the matrices $S$ in (84) and $\partial_{\mathbfit{R}}V$ in (133) and the block partitioning of the matrix $\partial_{\mathcal{M}}V$ in (134). $\blacksquare$

Note that the proof of Theorem 4 employs not only Lemma 3 but also the presence of the inverse matrix $\mathcal{P}^{-1}$ in (138) and (139) as a consequence of (132) due to the specific structure of the purity functional, leading to the self-cancellation of the matrix $S$ in

[TABLE]

Similarly to the steady-state covariances themselves, the relations (131) allow the infinitesimal perturbation analysis of the purity functional for a particular realization of the quantum system to be easily modified for an equivalent realization.

In view of the symplectic property of the matrices $S_{k}\in\mathrm{Sp}(\Theta_{k})$, the transformation (131) also shows that the partial Frechet derivatives $\rho_{k}$ and $\mu_{k}$ of the functional $V$ can not be made arbitrarily small. Moreover, since the norm of such matrices $S_{k}$ can be arbitrarily large, their inappropriate choice can make the purity functional highly sensitive to the parameters of the cascaded oscillators. In Section 8, a criterion will be specified for the minimization of this sensitivity, thus putting the choice of the symplectic matrices $S_{1},\ldots,S_{N}$ on a rational footing.

We will now apply Lemma 2 and elements of the proof of Theorem 4 to computing the logarithmic Frechet derivatives of the purity functional. The formulation of the following theorem employs the observability Gramian $\mathcal{Q}:=(\mathcal{Q}_{jk})_{1\leqslant j,k\leqslant N}$ of the pair $(\mathcal{A},\mathcal{P}^{-1/2})$, which is split into blocks $\mathcal{Q}_{jk}\in\mathbb{R}^{n_{j}\times n_{k}}$ and satisfies the ALE

[TABLE]

where, as before, $\mathcal{P}$ is the controllability Gramian of $(\mathcal{A},\mathcal{B})$ from (73). Also, use will be made of the Hankelian $\mathcal{H}:=(\mathcal{H}_{jk})_{1\leqslant j,k\leqslant N}$ for the triple $(\mathcal{A},\mathcal{B},\mathcal{P}^{-1/2})$, consisting of the blocks $\mathcal{H}_{jk}\in\mathbb{R}^{n_{j}\times n_{k}}$ and defined by

[TABLE]

The matrix $\mathcal{H}$ is diagonalizable and its eigenvalues are the squared Hankel singular values of the triple $(\mathcal{A},\mathcal{B},\mathcal{P}^{-1/2})$.

Theorem 5**.**

Suppose the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz. Then the Frechet derivatives (130) can be computed as

[TABLE]

for all $k=1,\ldots,N$. Here, $\mathcal{Q}_{\bullet k}:={\small\begin{bmatrix}\mathcal{Q}_{1k}\\ \vdots\\ \mathcal{Q}_{Nk}\end{bmatrix}}\in\mathbb{R}^{n\times n_{k}}$ denotes the $k$th block-column of the observability Gramian $\mathcal{Q}$ from (140), and $\mathcal{H}_{kk}$ is the $k$th diagonal block of the Hankelian $\mathcal{H}$ from (141). $\square$

Proof.

The representations of the partial Frechet derivatives in (108) and (109) allow their adjoint operators $\partial_{\mathbfit{R}}\mathcal{P}^{\dagger}$ and $\partial_{\mathcal{M}}\mathcal{P}^{\dagger}$ to be computed as

[TABLE]

where use is made of self-adjointness of the symmetrizer $\mathbf{S}$ (which is the orthogonal projection onto $\mathbb{S}_{n}$) together with the relation (137) which leads to $\mathbf{L}_{\alpha}^{\dagger}=\mathbf{L}_{\alpha^{\mathrm{T}}}$ and is combined with the symmetry of $\mathcal{P}$ and antisymmetry of $\mathbfit{\Theta}$. Here, $\mathbfit{\Sigma}$ denotes the orthogonal projection of $\mathbb{R}^{n\times n}$ onto the subspace $\mathfrak{R}$ of block-diagonal matrices $\mathbfit{R}$ in (69), which maps an arbitrary matrix $\zeta:=(\zeta_{jk})_{1\leqslant j,k\leqslant N}\in\mathbb{R}^{n\times n}$ with blocks $\zeta_{jk}\in\mathbb{R}^{n_{j}\times n_{k}}$ to the matrix

[TABLE]

Substitution of (144) and (146) into (133) leads to

[TABLE]

where use is made of the symmetric matrix $\mathcal{Q}=\mathbf{L}_{\mathcal{A}^{\mathrm{T}}}(\mathcal{P}^{-1})$ from (140) along with the Hankelian $\mathcal{H}$ from (141) and the block diagonal structure of the CCR matrix $\mathbfit{\Theta}$ in (65). The representation (142) is now obtained by considering the diagonal blocks of (147). By a similar reasoning, substitution of (145) into (134) leads to

[TABLE]

The structure of the operator $\Gamma$ in (106) implies that its adjoint $\Gamma^{\dagger}$ acts on a matrix $\alpha:=(\alpha_{jk})_{1\leqslant j,k\leqslant N}\in\mathbb{R}^{n\times n}$ with blocks $\alpha_{jk}\in\mathbb{R}^{n_{j}\times n_{k}}$ as

[TABLE]

where $\mathbf{T}$ denotes the operator of matrix transpose. Here, $\mathbfit{\Phi}$ is a linear operator which maps the matrix $\alpha$ to $\mathbfit{\Phi}(\alpha):=(\mathbfit{\Phi}_{jk}(\alpha))_{1\leqslant j,k\leqslant N}\in\mathbb{A}_{n}$ whose blocks are given by

[TABLE]

and depend only on the block lower triangular part of $\alpha$ (including its diagonal blocks). Also, use has been made of the following representation of the operators in (106):

[TABLE]

in combination with the antisymmetry of the matrix $J$ and the fact that both the antisymmetrizer $\mathbf{A}$ (which is an orthogonal projection) and $\mathbf{T}$ are self-adjoint operators.555More precisely, the transpose $\mathbf{T}_{r,s}$, applied to $(r\times s)$-matrices, satisfies $\mathbf{T}_{r,s}^{\dagger}=\mathbf{T}_{s,r}$. According to (150), the $k$th block of the image matrix in (149) is computed as

[TABLE]

for all $k=1,\ldots,N$. Therefore, by applying (151) to the matrix $\mathbfit{\Theta}\mathcal{H}$ in (148) and recalling the block diagonal structure and antisymmetry of $\mathbfit{\Theta}$ in (65), it follows that

[TABLE]

which establishes (143). $\blacksquare$

Due to the block lower triangular structure of the matrix $\mathcal{A}$ (that is, block upper triangularity of $\mathcal{A}^{\mathrm{T}}$), the ALE (140) of Theorem 5 for the observability Gramian $\mathcal{Q}$ can also be solved in a recursive fashion, similar to the computation of the controllability Gramian $\mathcal{P}$ in Lemma 1. However, in contrast to $\mathcal{P}$, the blocks of the matrix $\mathcal{Q}$ satisfy recurrence equations which unfold backwards (that is, from the last oscillator in the cascade towards the first one), thus being reminiscent of Bellman’s dynamic programming equations. Note that here the role of time is played by the spatial parameter which numbers the oscillators in the cascade.

7 Recursive computation of the Frechet derivatives

Although Theorem 5 provides a complete set of equations for computing the Frechet derivatives $\rho_{1},\ldots,\rho_{N}$ and $\mu_{1},\ldots,\mu_{N}$ of the purity functional in (130), we will outline its recursive version which takes into account the cascade structure of the system to a fuller extent. For what follows, the matrix pairs

[TABLE]

which specify the energetics of individual oscillators and take values in the corresponding spaces $\mathfrak{E}_{k}:=\mathbb{S}_{n_{k}}\times\mathbb{R}^{m\times n_{k}}$, are assembled into an $N$-tuple

[TABLE]

Accordingly, $E$ takes values in the Hilbert space $\mathfrak{E}:=\mathfrak{E}_{1}\times\ldots\times\mathfrak{E}_{N}$ with the direct-sum inner product $\langle\cdot,\cdot\rangle_{\mathfrak{E}}$ (generated from the Frobenius inner products of matrices). The Frechet derivative of $V$ (as a composite function of $E$) takes the form

[TABLE]

which combines (130), (133) and (134). Note that (154) is obtained without using the spatial causality mentioned in Section 5. This property can be taken into account in order to gain a computational advantage for long cascades. More precisely, the decomposition (97) of the functional $V$ and the fact that $V_{1},\ldots,V_{k-1}$ are independent of the matrix pair $E_{k}$ in (152) imply that

[TABLE]

for all $k=1,\ldots,N$, where

[TABLE]

is the “tail” part of $V$. Here, $\Pi_{\geqslant k}$ denotes the Schur complement of the block $\mathcal{P}_{k-1}$ in the matrix $\mathcal{P}$ in (87) computed as

[TABLE]

with

[TABLE]

according to the block partitioning

[TABLE]

which is similar to (88). The second equality in (157) provides a probabilistic interpretation of $\Pi_{\geqslant k}$ in terms of the auxiliary classical Gaussian random vectors $\xi_{1},\ldots,\xi_{N}$ of Section 5, with

[TABLE]

Accordingly, the blocks $Q_{\geqslant k}$ and $\mathcal{P}_{\geqslant k}$ in (159) are the covariance matrices for the random vectors $\Xi_{k-1}$ in (95) and $\Xi_{\geqslant k}$ in (160):

[TABLE]

Note that the first block-row of the matrix $Q_{\geqslant k}$ in (161) is the matrix $Q_{k}$ given by (89). Similarly to (99), the probabilistic meaning of the matrix $T_{k}$ in (158) is described in terms of classical conditional expectations as

[TABLE]

In particular, in view of (162),

[TABLE]

is the unconditional covariance matrix of the vector $\Xi_{\geqslant 1}=\Xi_{N}$. The subsequent Schur complements $\Pi_{\geqslant 2},\ldots,\Pi_{\geqslant N}=\Pi_{N}$ are amenable to a recursive computation as follows.

Lemma 6**.**

Suppose the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz. Then the Schur complement in (157) satisfies the recurrence equation

[TABLE]

for all $k=2,\ldots,N$, where the matrices

[TABLE]

are submatrices of the preceding Schur complement

[TABLE]

and the initial condition is given by (164). $\square$

Proof.

The relations (165)–(169) follow directly from the representation of the matrices in terms of the classical Gaussian random vectors $\xi_{1},\ldots,\xi_{N}$ of Section 5 and are similar to the conditional covariance matrix update in the discrete-time Kalman filter [1, 22]. $\blacksquare$

By a reasoning, similar to (132)–(134), the partial Frechet derivative $\partial_{E_{k}}V_{\geqslant k}$ of the tail functional $V_{\geqslant k}$ with respect to the matrix pair $E_{k}$ in (152) can be computed in terms of the Schur complement $\Pi_{\geqslant k}$ (see (155)–(157)) as

[TABLE]

For computing the Frechet derivative $\partial_{E_{k}}\Pi_{\geqslant k}$ in the following lemma, which is required for (170), we will use (159) along with the partitioning

[TABLE]

whose blocks are recovered from (58) and (59). More precisely,

[TABLE]

is a block lower triangular matrix with the diagonal blocks $A_{k},\ldots,A_{N}$, and

[TABLE]

where use is made of the matrices $\mathcal{C}_{k-1}$ from (60) and $\mathcal{B}_{\geqslant k}$ from (94). The following lemma, which is similar to Lemmas 1 and 2, carries out an infinitesimal perturbation analysis of the Schur complements in (157).

Lemma 7**.**

Suppose the matrices $A_{1},\ldots,A_{N}$ in (46) are Hurwitz. Then, for any $k=1,\ldots,N$, the Frechet derivatives of the Schur complement $\Pi_{\geqslant k}$ in (157) with respect to the energy and coupling matrices of the $k$th oscillator can be computed as

[TABLE]

Here, $\Pi_{k\bullet}\in\mathbb{R}^{n_{k}\times(n_{k}+\ldots+n_{N})}$ denotes the first block-row of the matrix $\Pi_{\geqslant k}$, and use is made of the operators (102)–(105) together with the matrix

[TABLE]

where $T_{k}$ is given by (158). $\square$

Proof.

Similarly to (90) of Lemma 1, for any $k=2,\ldots,N$, the block $Q_{\geqslant k}$ of the matrix $\mathcal{P}$ in (159) satisfies the ASE

[TABLE]

where the matrices $\mathcal{B}_{\geqslant k}$ and $D_{k}$ are given by (94) and (173). Similarly to (110) in the proof of Lemma 2, the first variation of the equation (177) with respect to $E_{k}$ takes the form

[TABLE]

Here, use is made of the spatial causality discussed above, whereby the matrices $\mathcal{A}_{k-1}$, $\mathcal{B}_{k-1}$, $\mathcal{P}_{k-1}$, pertaining to the first $k-1$ oscillators in the cascade, are independent of $E_{k}$, and hence, their variations vanish. Also, in view of (172), the first variations of the $E_{k}$-dependent matrices $\mathcal{A}_{\geqslant k}$, $\mathcal{B}_{\geqslant k}$, $D_{k}$ in (178) are computed as

[TABLE]

where the matrices $\mathcal{B}_{>k}$ and $\mathcal{C}_{k-1}$ are independent of $E_{k}$. In terms of the operators (102) and (104), the solution of the ASE (178) leads to

[TABLE]

Here, in view of (179)–(181), the components of the Frechet derivatives $\partial_{E_{k}}\mathcal{A}_{\geqslant k}$, $\partial_{E_{k}}\mathcal{B}_{\geqslant k}$ and $\partial_{E_{k}}D_{k}$ are given by

[TABLE]

Substitution of (183)–(188) into (182) leads to

[TABLE]

where use is made of the relation $\begin{bmatrix}I_{n_{k}}&0\end{bmatrix}Q_{\geqslant k}=Q_{k}$ which follows from (89) and (161). Since $\mathcal{P}_{k-1}$ does not depend on $E_{k}$, the corresponding derivatives of the matrix $T_{k}$ from (158) are expressed in terms of (189) and (190) as

[TABLE]

Now, in view of the partitioning (159), the matrix $\mathcal{P}$ can be factorized in terms of $T_{k}$ and the matrix $\Pi_{\geqslant k}$ from (157) as

[TABLE]

for any $k=1,\ldots,N$. This factorization is closely related to the decomposition

[TABLE]

where, in view of (163), the random vectors $\Xi_{k-1}$ and $\Xi_{\geqslant k}-\widehat{\Xi}_{\geqslant k}$ are uncorrelated and have covariance matrices $\mathcal{P}_{k-1}$ and $\Pi_{\geqslant k}$, respectively. Substitution of (192) into the ALE (73) allows the latter to be represented in the form

[TABLE]

where

[TABLE]

with the matrix $\widetilde{\mathcal{B}}_{\geqslant k}$ given by (176), and use is made of (171). By substituting (194) and (195) into (193), it follows that the second diagonal block of this ALE takes the form

[TABLE]

By a reasoning, similar to (182), the Frechet differentiation of both parts of (196) with respect to $E_{k}$ leads to

[TABLE]

Here,

[TABLE]

in view of (191), (176) and the fact that the matrix $\mathcal{B}_{k-1}$ does not depend on $E_{k}$. Substitution of (183)–(190) into (197) and (198) leads to

[TABLE]

which establishes (174) and (175) in view of the representation $\Pi_{k\bullet}=\begin{bmatrix}I_{n_{k}}&0\end{bmatrix}\Pi_{\geqslant k}$ for the first block-row of the matrix $\Pi_{\geqslant k}$. $\blacksquare$

Application of Lemma 7 to (170) is carried out similarly to the proof of Theorem 5 by using the relations

[TABLE]

Here, the adjoint operators can be found from (174) and (175) (the resulting expressions are cumbersome and are omitted). Their evaluation at $\Pi_{\geqslant k}^{-1}$ in (199) involves the observability Gramian $\mathcal{Q}_{\geqslant k}$ of the pair $(\mathcal{A}_{\geqslant k},\Pi_{\geqslant k}^{-1/2})$ satisfying the ALE

[TABLE]

which reduces to (140) in the case $k=1$. Note that the dimension of (200) decreases with $k$ and becomes relatively small at the end of the cascade $k=N$.

8 Minimization of the mean square sensitivity index

In what follows, we will use the vectorization of matrices [24, 56] in regard to the matrix-valued parameters of the cascaded oscillators and the corresponding Frechet derivatives. The full (rather than half-) vectorization is denoted by $\vec{(\cdot)}$ or $\mathrm{vec}(\cdot)$ interchangeably and is applicable to symmetric matrices and assemblages of matrices. Accordingly, the $N$-tuple $E$ of the pairs $E_{k}$ in (152) and (153) is vectorized as

[TABLE]

where $\vec{R}_{k}\in\mathbb{R}^{n_{k}^{2}}$, so that $\vec{E}_{k}\in\mathbb{R}^{n_{k}(n_{k}+m)}$. In terms of the vectorization, the linear operator in (102) is represented as

[TABLE]

where $\alpha\oplus\beta:=\alpha\otimes I+I\otimes\beta$ is the Kronecker sum of matrices. Furthermore, we denote by $\Upsilon_{r}$ the duplication $r^{2}\times\frac{1}{2}r(r+1)$-matrix [24, 56] which relates the full vectorization $\vec{M}$ of a symmetric matrix $M$ of order $r$ to its half-vectorization $\mathrm{vech}(M)$ (that is, the column-wise vectorization of its triangular part below and including the main diagonal) by

[TABLE]

For example,

[TABLE]

The presence of linear degeneracies in the full vectorization $\vec{E}_{k}$ in (201) (coming from the symmetry of the energy matrices $R_{k}\in\mathbb{S}_{n_{k}}$) can be taken into account by considering a smaller number of independent variables as

[TABLE]

Here, the matrix $\mho_{k}$ is of full column rank, and, in contrast to $\vec{E}_{k}$, the entries of the vector $e_{k}\in\mathbb{R}^{n_{k}(\frac{1}{2}(n_{k}+1)+m)}$ are linearly independent.

Now, suppose the energy and coupling matrices of the component oscillators are subject to infinitesimal perturbations (or implementation errors) $\delta R_{k}$ and $\delta M_{k}$, which, in accordance with (201) and (204), are represented in the vectorized form as $\delta\vec{E}_{k}$ or $\delta e_{k}$. The corresponding (linearized) variation of the functional $V$ is

[TABLE]

where

[TABLE]

and use is made of the vectorizations of the partial Frechet derivatives $\rho_{k}$ and $\mu_{k}$ from (130). If the modelling errors $\delta e_{1},\ldots,\delta e_{N}$ (which encode the vectors $\delta\vec{E}_{1},\ldots,\delta\vec{E}_{N}$) are of a classical random nature and are statistically uncorrelated with each other for different subsystems, they can be regarded as zero mean random vectors with covariance matrices

[TABLE]

for all $j,k=1,\ldots,N$. Here, $\epsilon>0$ is a small scaling parameter, and $\Sigma_{1},\ldots,\Sigma_{N}$ are appropriately dimensioned positive definite matrices (which specify the shape of the scattering ellipsoids). Note that (207) completely specifies the covariance operators $\mathbf{E}(\delta E_{j}\otimes\delta E_{k})$ for the random matrices $\delta E_{1},\ldots,\delta E_{N}$, whereby the variance of the quantity $\delta V$ in (205) takes the form

[TABLE]

and is proportional (with a constant coefficient $\epsilon$) to the functional

[TABLE]

Here, use is also made of the relations (155) and (156). Note that $Z$ in (209) is a quadratic function of the gradient $\partial_{\vec{E}}V$ and, in view of (208), can be regarded as a mean square measure for the sensitivity of the purity functional to the modelling errors in the energy and coupling matrices.

As mentioned in Section 6, the sensitivity index $Z$ is not invariant under the symplectic similarity transformations (83) of the component oscillators. Moreover, by applying (131) of Theorem 4, it follows that the gradient vectors in (206) are transformed as

[TABLE]

where the partial Frechet derivatives $\rho_{1},\ldots,\rho_{N}$ and $\mu_{1},\ldots,\mu_{N}$ are evaluated at the original realization of the system as described in Theorem 5 and Section 7. Substitution of (210) into (209) leads to the following dependence of the quantities $Z_{k}$ in (209) on the transformation matrices $S_{k}\in\mathrm{Sp}(\Theta_{k})$:

[TABLE]

Each of the functions $\Phi_{k}$ in (211) is a quartic polynomial of the entries of $S_{k}$, with the coefficients of this polynomial depending on the vector $d_{k}$ and the matrix $\Sigma_{k}$.

In the context of robust Gaussian state generation, the above discussion suggests that the symplectic transformation matrices $S_{1},\ldots,S_{N}$ (for a given cascade realization) can be found so as to minimize the mean square sensitivity index $Z$ in (209). Since $Z$ has an additive structure, its minimization can be split into $N$ independent problems of minimizing the functions $\Phi_{1},\ldots,\Phi_{N}$ in (211) over the corresponding symplectic groups:

[TABLE]

The (optimal or suboptimal) solutions $S_{1},\ldots,S_{N}$ of these problems can be used for balancing the original realization of the cascaded oscillators. The resulting equivalent realization of the composite quantum system will have a decreased mean square sensitivity of the purity functional with respect to small random perturbations. This criterion differs from the previously known balanced realizations of classical linear systems and their quantum analogues [33] which equate the controllability and observability Gramians of the system (in the spirit of the Kalman duality principle). Another related (entropy theoretic) optimality criterion will be discussed in Appendix A.

We will now assume that there is prior information on parametric uncertainties which ensures the following upper bounds on the covariance matrices $\Sigma_{k}$ in (207):

[TABLE]

where $a_{k}$ and $b_{k}$ are positive scalars for all $k=1,\ldots,N$. In particular, (213) holds if the uncertainties $\delta R_{k}$ and $\delta M_{k}$ in the energy and coupling matrices are uncorrelated and satisfy

[TABLE]

with $\epsilon$ the small scale factor as before. Since the operator norm of the duplication matrix $\Upsilon_{r}$ in (203) does not exceed $\sqrt{2}$ for any $r$, then (213) implies that

[TABLE]

where the matrix $\mho_{k}$ is given by (204). The last inequality in (215) leads to the following upper bound on the function $\Phi_{k}$ in (211):

[TABLE]

where we have used the isometric property $|\vec{M}|=\|M\|$ of the full vectorization for any matrix $M$. Therefore, a suboptimal solution to the problem of minimizing the function $\Phi_{k}$ in (212) can be obtained by replacing it with the right-hand side of (216):

[TABLE]

This suboptimal approach can be applied to all the oscillators in the cascade by solving (217) independently for every $k=1,\ldots,N$. In the case $n_{k}=2$, when $O_{k}$ is a one-mode oscillator and the CCR matrix $\Theta_{k}$ in (43) coincides with the matrix $\mathbf{J}$ from (7) up to a scalar factor, the replacement problem reduces effectively to a quadratic optimization problem and its solution is provided below.

Theorem 8**.**

Suppose $\rho\in\mathbb{S}_{2}$ and $\mu\in\mathbb{R}^{m\times 2}$ are given matrices, the second of which is of full column rank:

[TABLE]

Then any locally optimal solution of the minimization problem

[TABLE]

can be represented as

[TABLE]

Here,

[TABLE]

is the matrix of rotation by an arbitrary angle $\phi$, and $U\in\mathbb{S}_{2}$ is a positive definite matrix which is computed as

[TABLE]

where the function

[TABLE]

of a scalar variable $z$ is applied to the real symmetric matrix $\tau^{-1/2}\rho\tau^{-1/2}$ and depends on a parameter $\lambda>0$. The latter is found as a unique solution of the equation

[TABLE]

where $r_{1}$ and $r_{2}$ are the eigenvalues of $\rho\tau^{-1}$. $\square$

Proof.

In view of the symmetry of the matrices $\rho$ and $\tau$ in (218), the quartic polynomial $\Psi$ in (219) can be represented as

[TABLE]

which is a quadratic function of a real positive semi-definite symmetric matrix

[TABLE]

Due to the identity $\sigma\mathbf{J}\sigma^{\mathrm{T}}=\det\sigma\mathbf{J}$, which holds for any matrix $\sigma\in\mathbb{R}^{2\times 2}$, the symplectic property $\sigma\in\mathrm{Sp}(\mathbf{J})$ is equivalent to $\det\sigma=1$, and hence, the problem (219) reduces to a constrained quadratic minimization problem

[TABLE]

By endowing the constraint $\ln\det U=0$ with a Lagrange multiplier $\lambda\in\mathbb{R}$, the Lagrange function for (227) takes the form

[TABLE]

where the $\frac{1}{2}$ factor is introduced for convenience. Any locally optimal solution $U$ of the problem (227) is necessarily a stationary point of the Lagrange function $\mathcal{L}_{\lambda}(U)$ for some $\lambda$, and hence, the corresponding (unconstrained) Frechet derivative of (228) vanishes:

[TABLE]

where use is made of (225) and (132). Since the matrices $\tau$ and $U$ in (218) and (226) are positive definite, and $\rho$ is symmetric, then (229) implies that $\lambda>0$ and

[TABLE]

By left and right multiplying both sides of this equality by $\sqrt{\tau}$, it follows that an appropriate transformation of the matrix $U$ satisfies

[TABLE]

where

[TABLE]

is an auxiliary real symmetric matrix. A reasoning, similar to that for solving a class of algebraic Riccati equations in [64, Lemma 10.1], shows that (230) has a unique solution $T\succ 0$ and this solution is representable as

[TABLE]

Here, for any given $\lambda>0$, a function $f_{\lambda}(z)$ of a real variable $z$ is evaluated at the real symmetric matrix $\widetilde{\rho}$ (see, for example, [13]), whereby the matrix $T$ in (232) commutes with $\widetilde{\rho}$. This allows $f_{\lambda}(z)$ to be found as the unique positive solution of the scalar equation

[TABLE]

or, equivalently, as the unique positive root $w$ of the quadratic polynomial $z^{2}w^{2}+w-\frac{\lambda}{2}$, which leads to (223). Now, since $\det\sigma=1$, then (226), (230) and (232) imply that the Lagrange multiplier $\lambda$ must satisfy

[TABLE]

where $r_{1}$ and $r_{2}$ are the eigenvalues of the matrix $\widetilde{\rho}$ in (231). The right-hand side of (233) is a strictly increasing continuous function of $\lambda>0$, which varies from 0 to $+\infty$, whereby (224) has a unique solution $\lambda$. The corresponding matrix $U$ in (230) is recovered from $T$ in (232) as described by (222), and the symplectic matrix $\sigma$ is found from (226) according to (220) and (221). It now remains to note that the matrix $\widetilde{\rho}=\tau^{-1/2}\rho\tau^{-1}\sqrt{\tau}$ in (231) is obtained through a similarity transformation from and is, therefore, isospectral to $\rho\tau^{-1}$. $\blacksquare$

In addition to its monotonicity, the left-hand side $h(\lambda)$ of (224) is a smooth convex function of $\lambda$ (see, for example, Fig. 3), which makes this equation amenable to effective numerical solution by a standard root-finding algorithm. In particular, application of the Newton-Raphson method with the initial condition

[TABLE]

produces a monotonically decreasing sequence $\lambda_{1},\lambda_{2},\ldots$ as $\lambda_{k+1}:=\lambda_{k}+\frac{\det\tau-h(\lambda_{k})}{h^{\prime}(\lambda_{k})}$ (so that $\lambda_{0}\leqslant\lambda_{1}\geqslant\lambda_{2}\geqslant\ldots$) which converges at a quadratic rate to the solution of (224), with the derivative of the function $h$ computed as

[TABLE]

in view of (223). The particular choice of the initial condition in (234) yields $h(\lambda_{0})\leqslant\frac{1}{4}\lambda_{0}^{2}=\det\tau$, whereby $\lambda_{0}$ is a guaranteed lower bound for the solution.

Although the “quartic-to-quadratic” reduction (225) remains valid for many-mode settings, the proof of Theorem 8 relies on the surjectivity of the map $\sigma\mapsto U$ in (226) in the sense that any positive definite matrix $U\in\mathbb{S}_{2}$ with $\det U=1$ is factorizable as $U=\sigma^{\mathrm{T}}\sigma$ for some symplectic matrix $\sigma\in\mathrm{Sp}(\mathbf{J})$. This property does not extend to higher dimensions $\nu$, for which the image of the corresponding symplectic group is a proper subset of $\mathfrak{U}:=\{U\in\mathbb{S}_{\nu}:\ U\succ 0,\ \det U=1\}$. In this case, the minimization of the quadratic function $\widetilde{\Psi}(U)$ in (227) over the larger set $\mathfrak{U}$ can only provide a lower (and not necessarily achievable) bound for the actual minimum value of the quartic polynomial $\Psi(\sigma)$.

9 An application to balancing cascaded one-mode oscillators

Consider a cascade of one-mode OQHOs with dimensions

[TABLE]

In regard to such systems, it can be assumed, without loss of generality, that the CCR matrices in (43) are given by

[TABLE]

In accordance with (20) and (7), this corresponds to the pairs of position and momentum operators as system variables of the component oscillators. The symplectic similarity transformations of the oscillators are specified by elements $S_{1},\ldots,S_{N}$ of the common symplectic group $\mathrm{Sp}(\mathbf{J})$ (regardless of the factor $\frac{1}{2}$ in (236)) which can be represented as

[TABLE]

in terms of positive scalars $\varsigma_{k}$ and angles $\phi_{k}$, $\psi_{k}$ (see, for example, [68]), where use is made of the rotation matrix (221).

In the framework of the mean square sensitivity minimization for the purity functional, described in Section 8, the matrix $S_{k}$ can be found by independently solving the optimization problem (217), associated with the $k$th oscillator, for every $k=1,\ldots,N$. In the one-mode case (235), this amounts to applying Theorem 8 with the parameters

[TABLE]

where $\rho_{k}$ and $\mu_{k}$ are the Frechet derivatives in (130) evaluated at the original realization of the component oscillator, and $a_{k}$ and $b_{k}$ are the upper bounds for the covariances of the uncertainties in the energy and coupling matrices in the sense of (213) (or (214) if they are uncorrelated). A comparison of (237) with (220) leads to

[TABLE]

Therefore, whilst the angle $\phi_{k}$ in (237) can be arbitrary, $\varsigma_{k}$ and $\frac{1}{\varsigma_{k}}$ are the eigenvalues of the corresponding matrix $U$, specified in Theorem 8, and the other angle $\psi_{k}$ is such that the columns of the matrix $\mathsf{R}(\psi_{k})^{\mathrm{T}}=\mathsf{R}(-\psi_{k})$ are the eigenvectors of $U$.

We will now provide the results of a numerical experiment on balancing a cascade of $N=3$ one-mode oscillators with $m=6$ input field channels and the energy and coupling matrices

[TABLE]

which were randomly generated subject to the condition that the matrices $A_{1}$, $A_{2}$, $A_{3}$ in (46) are Hurwitz. The Frechet derivatives

[TABLE]

in (130) are computed as described in Section 7. The upper bounds $a_{k}$ and $b_{k}$ (also randomly generated) on the covariances of the parametric uncertainties in (213) and the contributions $\Psi_{k}(S_{k})$ from the component oscillators in (216) to the mean square sensitivity index before and after the balancing (82) (and their ratios) are given in Tab. 1.

The symplectic transformation matrices, which minimize the functions $\Psi_{k}$ for the oscillators, are found by using Theorem 8 in combination with (238), (239) and the Newton-Raphson algorithm specified in Section 8:

[TABLE]

The reduction in the total value of the mean square sensitivity index, which is achieved for these matrices, is

[TABLE]

which is an intermediate value in comparison with the reduction ratios for the component oscillators at the bottom of Tab. 1.

10 Conclusion

We have considered an approach to balancing cascaded linear quantum stochastic systems via symplectic similarity transformations of their realizations for the generation of Gaussian invariant quantum states. The optimality criterion considered is based on minimizing the mean square sensitivity of the purity functional quantified by the norm of its logarithmic Frechet derivative with respect to the energy and coupling matrices subject to small random perturbations. We have discussed a recursive computation of this performance functional together with a link to the classical Fisher information metric and reduced its optimization to a set of independent problems of minimizing quartic polynomials on symplectic groups. We have developed an effective algorithm for solving these problems in the one-mode setting when they reduce to quadratic optimization problems. We have also outlined the infinitesimal perturbation analysis for translation invariant cascades of identical oscillators using spatial $z$-transforms. The results of this paper may be of use in the robust generation of pure Gaussian states for the purposes of state preparation in quantum computing and quantum information processing involving quantum stochastic networks, especially those with fractal-like architectures (in the form of acyclic directed graphs such as trees). Beyond the quantum domain, similar ideas are applicable to perturbation analysis and optimization of stochastic versions of port-Hamiltonian systems [59] whose dynamics bear special structures reflecting the energetics of classical physical systems.

Appendix A Classical relative entropy deviation for Gaussian states

We will now discuss a connection of the mean square sensitivity index (209) for the purity functional with an alternative criterion which quantifies the perturbations in the matrix $\mathcal{P}$ in terms of classical distance measures applied to invariant Gaussian quantum states.666As mentioned in Section 4, inaccuracies in the energy and coupling matrices do not affect the CCRs (66), thus allowing the Gaussian invariant states to be distinguished by the real parts of their quantum covariance matrices. In particular, an important role is played by entropy theoretic statistical distances. Recall that the relative entropy (or the Kullback-Leibler informational divergence [4]) of a probability measure $\mathsf{P}$ with respect to a reference probability measure $\mathsf{P}_{*}$ (on the same measurable space $(\mathfrak{X},\mathfrak{B})$) is defined as the expectation

[TABLE]

Here, $\mathsf{P}$ is assumed to be absolutely continuous with respect to $\mathsf{P}_{*}$, and $\frac{\mathrm{d}\mathsf{P}}{\mathrm{d}\mathsf{P}_{*}}$ denotes the corresponding Radon-Nikodym derivative [50, 51]. The last equality in (240) links the functional $\mathbf{D}$ with the relative Renyi $\alpha$-entropy [46] which, in the case $\alpha=\frac{1}{2}$, is present in the Hellinger distance between $\mathsf{P}$ and $\mathsf{P}_{*}$ given by

[TABLE]

While the choice of a statistical distance is not critical for our purposes, the Kullback-Leibler relative entropy (240) is particularly convenient in application to Gaussian distributions and has direct links with the Shannon information theory [12].

The following lemma, which is concerned with the Gaussian case, is well-known in the context of the maximum entropy principle [4, 25] (see also [41, Lemma 4 on pp. 313–314], [64, Lemma 9.1 on p. 1129] and references therein).

Lemma 9**.**

For two zero-mean Gaussian distributions $\mathcal{N}(0,P)$ and $\mathcal{N}(0,P_{*})$ in $\mathbb{R}^{n}$ with nonsingular covariance matrices $P$ and $P_{*}$, the relative entropy (240) of $\mathcal{N}(0,P)$ with respect to $\mathcal{N}(0,P_{*})$ takes the form

[TABLE]

where

[TABLE]

The asymptotic behaviour of this quantity, as $P\to P_{*}$, is described by

[TABLE]

$\square$**

The relation (243)777which follows from (241) and (242) in view of the Frechet derivatives (132) and $\partial_{\chi}^{2}\ln\det\chi(M)=-\chi^{-1}M\chi^{-1}$ on the set of real positive definite symmetric matrices $\chi$. shows that, in a small neighbourhood of the reference matrix $P_{*}$, the quantity $\sqrt{\mathbfit{\Delta}(P,P_{*})}$ gives rise to the relative deviation of the actual covariance matrix $P$ from $P_{*}$:

[TABLE]

Here, $g$ is a positive definite self-adjoint operator which acts on the Hilbert space $\mathbb{S}_{n}$ as

[TABLE]

and specifies the local metric tensor associated with the Fisher information distance [48] (on the set of positive definite covariance matrices of order $n$ regarded as a Riemannian manifold).

In application to the robustness of Gaussian state generation with respect to the implementation errors (modelled in Section 8 as zero-mean classical random vectors with the covariance matrices (207)), the relations (244) and (245) suggest that the balancing of the oscillators can also be based on minimizing the functional

[TABLE]

where $\mho_{k}$ is the matrix from (204). Indeed, averaging of the leading quadratic term of the Kullback-Leibler deviation of the perturbed covariance matrix from its nominal value leads to

[TABLE]

which is proportional to the quantity $Z$ in (246) with a constant multiplier $\epsilon$. In (247), use has also been made of the symmetry of the matrix $\mathcal{P}$ together with the identity $\mathcal{P}^{-1}\otimes\mathcal{P}^{-1}=(\mathcal{P}\otimes\mathcal{P})^{-1}$. Note that the entropy theoretic criterion (246) yields an upper bound on the purity functional sensitivity index (209) in the sense that

[TABLE]

in view of the Cauchy-Bunyakovsky-Schwarz inequality (applied here to the Frobenius inner product), though the coefficient of proportionality $n$ in (248) is large for many-mode systems. The above described connection is a consequence of the linear relation between the differential entropy [4] of a classical Gaussian distribution and the log determinant of its covariance matrix. In fact, the quantities $V_{k}$ and $V_{\geqslant k}$ in (97) and (156) (which are parts of the purity functional) are linearly related to the conditional differential entropies [4, 12] of the auxiliary classical Gaussian random vectors $\xi_{k}$ and $\Xi_{\geqslant k}$ in (160) with respect to the vector $\Xi_{k-1}$ in (95).

Appendix B Translation invariant cascades

We will now briefly discuss a translation invariant case when all the component oscillators, described in Section 3, are identical and have a common dimension $n$, CCR matrices $\Theta\in\mathbb{A}_{n}$, energy matrices $R\in\mathbb{S}_{n}$, coupling matrices $M\in\mathbb{R}^{m\times n}$, and the state-space matrices

[TABLE]

corresponding to (46)–(48), and the transfer functions in (76) and (77) given by

[TABLE]

In this case, the initial spaces of the oscillators are copies of a common Hilbert space $\mathfrak{H}$, and the QSDEs (44) and (45) take the form

[TABLE]

with $Y_{0}:=W$ as before. Following [63] (see also [53, 54]) and regarding the composite system as an infinite cascade of oscillators, we will use the spatial $z$-transforms of the system and output variables:

[TABLE]

These series are mean square convergent for sufficiently large values of the complex parameter $z$. More precisely, since the output variables of the oscillators are linearly related to the system variables and the input fields, the convergence of both series is guaranteed if

[TABLE]

provided the upper limit is finite, which is obtained by using an appropriate modification of the Cauchy-Hadamard theorem. The vectors $\mathcal{X}_{z}$ and $\mathcal{Y}_{z}$ consist of time-varying (and not necessarily self-adjoint) operators on the system-field space $\mathfrak{H}^{\otimes\infty}\otimes\mathfrak{F}$ and satisfy the CCRs

[TABLE]

which follow from the bilinearity of the commutator and the CCRs (43), provided $|z||v|>1$.

Lemma 10**.**

For all $z\in\mathbb{C}$ large enough (for example, such as in (257)), the processes $\mathcal{X}_{z}$ and $\mathcal{Y}_{z}$ in (255) and (256) satisfy the QSDEs

[TABLE]

Here, the $z$-dependent state-space matrices are computed as

[TABLE]

$\square$**

Proof.

Application of the $z$-transforms (255) and (256) to the sequence of QSDEs (253) and (254) with $k=1,2,\ldots$ leads to

[TABLE]

Here, use is also made of the convention $Y_{0}=W$ and the identities

[TABLE]

The QSDE (260) with the matrices (262) is obtained by solving (264) for $\mathrm{d}\mathcal{Y}_{z}$ as

[TABLE]

and substituting the result into (263), which leads to the QSDE (259) with the matrices (261). $\blacksquare$

In application to the matrices $A$ and $B$ of an individual oscillator in (249) and (250), the PR property (51) takes the form

[TABLE]

and leads to a similar PR condition for the matrices $\mathcal{A}_{z}$ and $\mathcal{B}_{z}$ in (261):

[TABLE]

which is equivalent to the preservation of the CCRs (258) under the QSDE (259).

Now, if the matrix $A$ in (249) is Hurwitz, then so is the matrix $\mathcal{A}_{z}$ in (261) for all sufficiently large $z$, and the set

[TABLE]

is nonempty and open. In this case, not only every finite subsystem of the infinite cascade has a unique Gaussian invariant state, but a similar property also holds for the process $\mathcal{X}_{z}$ which involves all the system variables of the cascade.

In view of (75), (251) and (252), the transfer function from the input quantum Wiener process $W$ to the vector $X_{k}$ of system variables of the $k$th oscillator in the translation invariant cascade is $FG^{k-1}$ for all $k=1,2,\ldots$. Hence (see also (78) and (79)), the transfer function from $W$ to $\mathcal{X}_{z}$ can be computed in the right-half plane $\mathrm{Re}s>0$ in terms of the transfer functions $F$ and $G$ and directly from the QSDE (259) as

[TABLE]

for any $z\in\mathbb{C}$ satisfying

[TABLE]

Therefore, if $\mathcal{A}_{z}$ is Hurwitz (that is, $z\in\mathfrak{Z}$ in view of (266)), the last equality in (267) implies that $\Phi_{z}$ belongs to the Hardy space $\mathcal{H}_{2}$. On the other hand, under the assumption that the matrix $A$ is Hurwitz, the second to last equality in (267) leads to an upper bound for the $\mathcal{H}_{2}$-norm of $\Phi_{z}$:

[TABLE]

in view of (268). In (269), we have also used the inequality888well known in the context of Banach algebras and their applications to the small gain theorem $\|(I_{m}-H)^{-1}\|_{\infty}\leqslant\frac{1}{1-\|H\|_{\infty}}$ which holds under the condition $\|H\|_{\infty}<1$ and follows from the submultiplicativity of the $\mathcal{H}_{\infty}$-norm. In particular, if, for given $z$ and $v$ from the set (266), the corresponding processes $\mathcal{X}_{z}$ and $\mathcal{X}_{v}$ have finite initial second moments (for example, both $z$ and $v$ satisfy (257) at time $t=0$), then they also have a steady-state quantum covariance matrix

[TABLE]

which is the unique solution of the ASE

[TABLE]

Note that (265) is part of this equation. In accordance with (87), the matrix $P_{jk}$ in (270) is the real part of the invariant quantum covariance matrix of the vectors $X_{j}$ and $X_{k}$ of system variables for the $j$th and $k$th oscillators. In accordance with (202), the vectorization of the ASE (271) allows its solution to be represented as

[TABLE]

In view of the rational dependence of the matrices $\mathcal{A}_{z}$ and $\mathcal{B}_{z}$ on $z$ in (261), it follows from (272) that $\mathcal{P}_{zv}$ is also a rational function of $z$ and $v$:

[TABLE]

Here, the matrix $A\oplus A$ is nonsingular (moreover, Hurwitz) due to $A$ being Hurwitz, and use is made of the auxiliary matrices

[TABLE]

The representation (273)–(275) can, in principle, be employed as a generating function of the steady-state matrices $P_{jk}$ for the infinite cascade. However, its straightforward use is complicated by the fact that, in general, the matrices $K$ and $L$ do not commute. We will therefore consider an upper bound for the covariance matrices.

Theorem 11**.**

Suppose the matrix $A$ of a component oscillator in the translation invariant cascade is Hurwitz. Then, for any $k=1,2,\ldots$, the steady-state covariance matrix of the system variables in the $k$th OQHO satisfies

[TABLE]

where $\|F\|_{2}$ and $\|G\|_{\infty}$ are the $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$-norms of the transfer functions in (251) and (252). $\square$

Proof.

A translation invariant version of the frequency-domain representation (74)–(77) for the steady-state covariances implies that

[TABLE]

for all $j,k=1,2,\ldots$. By using the spectral radius $\lambda_{\max}(\Omega)=2$ of the quantum Ito matrix $\Omega$ from (49) and submultiplicativity of the $\mathcal{H}_{\infty}$-norm, it follows from (277) that

[TABLE]

The inequality (276) can now be obtained by taking the trace on both sides of (278) in view of the antisymmetry of the CCR matrix $\Theta$. $\blacksquare$

In accordance with (80), the $\mathcal{H}_{\infty}$-norm of the transfer function $G$ of an individual oscillator in (252) satisfies $\|G\|_{\infty}\geqslant 1$. Therefore, (276) can only bound the exponential growth of the steady-state covariances for distant oscillators along the cascade. Also note that, since the operator norms of the matrices $P_{jj}$, $P_{jk}$, $P_{kk}$ satisfy

[TABLE]

(due to the positive semi-definiteness which the matrix $\small{\begin{bmatrix}P_{jj}&P_{jk}\\ P_{kj}&P_{kk}\end{bmatrix}}$ inherits from $\mathcal{P}_{N}$ in (87)), the upper bound (276) ensures absolute convergence of the series in (270) for all $z$ and $v$ satisfying (268), that is, $|z|,|v|>\|G\|_{\infty}$.