A Transverse Hamiltonian Approach to Infinitesimal Perturbation Analysis of Quantum Stochastic Systems

This paper is concerned with variational methods for open quantum systems with Markovian dynamics governed by Hudson–Parthasarathy quantum stochastic differential equations. These QSDEs are driven by quantum Wiener processes of the external bosonic fields and are specified by the system Hamiltonian and system–field coupling operators. We consider the system response to perturbations in these operators and introduce a transverse Hamiltonian which encodes the propagation of the perturbations through the unitary system–field evolution. This approach provides an infinitesimal perturbation analysis tool which can be used for the development of optimality conditions in quantum control and filtering problems. As performance criteria, such settings employ quadratic (or more complicated) cost functionals of the system and field variables to be minimized over the energy and coupling parameters of system interconnections. We demonstrate an application of the transverse Hamiltonian variational approach to a mean square optimal coherent quantum filtering problem for a measurement-free field-mediated cascade connection of a quantum system with a quantum observer.


Introduction
In comparison with classical mechanics which considers macroscopic objects both in deterministic and stochastic settings, quantum mechanics is concerned with physical phenomena at atomic and subatomic scales and inherently incorporates randomness. In particular, the squared absolute value of the wave function of a quantum mechanical particle is interpreted as a probability density function, mixed quantum states are built of pure ones using randomization, and the latter is present in the model of quantum measurement [1,2]. However, in contrast to scalar-valued classical probability measures [3], quantum probability describes statistical properties of quantum variables by using quantum states in the form of density operators [4,5] on the same Hilbert space where those variables act as linear operators.
The noncommutativity and canonical commutation structures originating from the operator-valued nature of quantum variables and quantum states give rise to specific features of quantum probability such as the absence of a classical joint probability distribution and conditional expectations for a set of noncommuting quantum variables (whereas an individual self-adjoint operator has a well-defined marginal distribution). Furthermore, because the microscopic realm is less amenable to manipulation by conventional macroscopic tools (unlike, for example, coin tossing as a manageable "random number generator" for thought and practical experiments in classical probability theory), its natural time evolution makes the statistical properties of quantum systems particularly tied to their dynamics. While an isolated quantum system undergoes a reversible evolution according to a one-parameter unitary group generated by the system Hamiltonian, a more realistic open dynamics scenario [6] involves interaction of the system with its environment, which can include other quantum or classical systems, measuring devices, and external quantum fields such as quantized electromagnetic radiation.
The interaction of an open quantum system with its surroundings is accompanied by energy exchange between the subsystems, and gives rise to dissipative effects. This provides a way to influence statistical and dynamic properties of such systems by arranging them into quantum networks [7] and varying the energy and coupling parameters of the resulting interconnections, which can consist of a quantum plant coupled in a measurement-based or coherent (that is, measurement-free) fashion to a controller or observer. This paradigm is used in quantum control [8][9][10], which develops systematic methods for achieving stability, robustness with respect to unmodelled dynamics, and optimality in the sense of relevant performance criteria for quantum systems and their applications such as quantum optics [11,12] and quantum information processing [13].
An important part in these developments belongs to quantum filtering, which, similarly to its classical predecessor (see for example [14,15]), is concerned with mean square optimal real-time estimation of the internal variables of a quantum plant using a measurementbased quantum Kalman filter [16][17][18] or a qualitatively different coherent quantum observer [19,20]. The latter does not process classical observations and, in contrast to the classical filters, does not compute the conditional expectations of the plant variables (this conditioning and related Bayesian approaches are not applicable in the noncommutative quantum case, as mentioned above). Instead, the coherent quantum filter is driven by the output quantum fields of the plant to produce a quantum process, and its performance as an estimator of the plant variables can be optimized in the sense of minimizing the mean square value of the "estimation error" by varying the parameters of the Hamiltonian and coupling operators of the filter.
The optimization problems arising in coherent quantum filtering and control (with the latter considering more complicated feedback interconnections of a quantum plant and a quantum controller) involve physical realizability (PR) constraints [21,22] which originate from the canonical commutation structures of quantum dynamic variables, the energetics of open quantum systems, and unitarity in the augmented system-environment evolution.
These features of open quantum dynamics are incorporated in Hudson-Parthasarathy quantum stochastic calculus [23,24] (see also [25]), which is often employed as a unified modelling framework in quantum filtering and control. In this approach, the external fields at the input of a quantum system of interest are represented by a multichannel quantum Wiener process with noncommuting components, which act on a symmetric Fock space and drive the quantum stochastic differential equation (QSDE) governing the system dynamics. In contrast to the classical SDEs with a standard Wiener process [26], the QSDE reflects the unitarity of the augmented system-field evolution, and its drift and dispersion terms (as well as the generator of Markovian quantum dynamics) are specified by the Hamiltonian and coupling operators. These operators, together with a scattering matrix which represents the photon exchange between the fields [24], describe the energetics of the quantum system and its interaction with the environment, and are usually functions (for example, polynomials or Weyl quantization integrals [27]) of the system variables. A particular form of this dependence and the commutation structure of the system variables affect the tractability of the quantum system under consideration.
In particular, quadratic dependence of the system Hamiltonian and linear dependence of the system-field coupling operators on the quantum mechanical position-momentum variables [28] lead to linear QSDEs for open quantum harmonic oscillators (OQHOs) [11,18], which play the role of building blocks in linear quantum control theory [10,21,[29][30][31][32]. The dynamics of such systems are relatively well understood and are similar to the classical linear SDEs in a number of respects, including the preservation of the Gaussian nature of quantum states [33,34] in the case of vacuum input fields. However, the coherent quantum analogue [35] of the classical linear-quadratic Gaussian (LQG) control problem [36,37] for OQHOs is complicated by the above mentioned PR constraints on the state-space matrices of the quantum controller and by the impossibility of taking advantage of the classical estimation-actuation separation principle and classical conditional expectations with their variational properties [15].
One of the existing approaches to the coherent quantum LQG (CQLQG) control and coherent quantum filtering (CQF) problems employs the Frechet derivatives of the mean square cost (with respect to the state-space matrices subject to the PR constraints) in obtaining the optimality conditions [20,38] and for numerical optimization [39]. This approach takes into account the quantum nature of the underlying problem only through the PR constraints, being "classical" in all other respects, which has advantages from the viewpoint of relying on well-developed conventional optimization methods. However, a disadvantage of this approach is that it is limited to certain parametric classes of linear controllers and observers. In particular, the resulting optimality conditions do not provide information on whether nonlinear quantum controllers or observers can outperform linear ones for linear quantum plants. For this reason, the coherent quantum control and filtering problems require novel variational methods for their solution, which would be able to operate with sensitivity of the system dynamics and relevant cost functionals to perturbations over wider classes of the Hamiltonian and coupling operators in a "coordinate-free" fashion.
To this end, the present paper (several of its results were briefly announced in [40]) outlines a fully quantum variational method which allows the sensitivity of the internal and output variables of a nonlinear quantum stochastic system to be investigated with respect to arbitrary (that is, not only linear-quadratic) perturbations of the Hamiltonian and coupling operators. This approach is based on using a transverse Hamiltonian, defined as an auxiliary time-varying self-adjoint operator which encodes the propagation of such perturbations through the unitary system-field evolution. More precisely, the perturbation of the quantum stochastic flow, which describes the time evolution of a system operator (a function of the system variables), is expressed in terms of the commutator with the transverse Hamiltonian. The resulting derivative processes for system operators lead to an infinitesimal perturbation formula for quantum averaged performance criteria (such as the mean square cost functional) which is applicable to the development of optimality conditions in coherent quantum control and filtering problems over larger classes of controllers and observers. In particular, this approach allows the sensitivity of OQHOs with quadratic performance criteria to be studied with respect to general perturbations of the Hamiltonian and coupling operators. We demonstrate its application to a CQF problem for a field-mediated cascade connection of a quantum plant with a quantum observer. In fact, the transverse Hamiltonian method has already been employed in [41] to establish the local sufficiency of linear observers in the mean square optimal CQF problem [20] for linear quantum systems with respect to varying the Hamiltonian and coupling operators of the observer along linear combinations of the Weyl operators [27]. As an extended version of the conference paper [40], the present work provides its results along with detailed proofs. Furthermore, as part of the additional material we discuss the joint commutation structure, including the cross-commutation relations, of the unperturbed system variables and their perturbations (see Section 3), and elaborate on the case where the perturbations of the Hamiltonian and coupling operators are in the Weyl quantization form (see Section 6). Note that our approach is different from [42], which develops a quantum Hamilton-Jacobi-Bellman principle for the density operator instead of the dynamic variables in a measurement-based quantum feedback control problem. We also note a parallel between the perturbation analysis discussed in the present paper and the fluctuation-dissipation theorem [43].
The rest of this paper is organized as follows. Section 2 specifies the class of quantum stochastic systems under consideration. Section 3 discusses the sensitivity of the internal and output variables of OQHOs to parametric perturbations within the families of quadratic system Hamiltonians and linear system-field coupling operators. Section 4 returns to general QSDEs and introduces the transverse Hamiltonian associated with arbitrary perturbations of the Hamiltonian and coupling operators of the system. The transverse Hamiltonian is used in Section 5 for infinitesimal perturbation analysis of system operators. Section 6 extends the transverse Hamiltonian method to quantum averaged performance criteria. Section 7 applies this approach to a mean square optimal CQF problem. Section 8 makes concluding remarks.

Quantum Stochastic Systems Being Considered
We consider an open quantum system which interacts with an external multichannel bosonic field and is equipped with dynamic variables X 1 (t), . . . , X n (t) evolving in time t 0 and assembled into a vector X(t) := (X k (t)) 1 k n (vectors are organized as columns unless indicated otherwise). These system variables are assumed to be self-adjoint operators on a composite system-field Hilbert space H ⊗ F , where ⊗ is the tensor product of spaces or operators (including the Kronecker product of matrices). Here, H is the initial space of the system, which provides a domain for X 1 (0), . . . , X n (0), and F is a symmetric Fock space [24] for the action of an even number m of quantum Wiener processes W 1 (t), . . . , W m (t). The latter are time-varying self-adjoint operators, which model the external fields and are assembled into a vector W(t) := (W k (t)) 1 k m . Unlike the classical Brownian motion [26] in R m , the quantum Wiener process W consists of noncommuting operator-valued components, and has a complex positive semi-definite Hermitian Ito matrix Ω := (ω jk ) 1 j,k m (identified with its tensor product Ω ⊗ I F with the identity operator I F on the Fock space F ) for its future-pointing Ito increments dW: Here, the transpose (·) T acts on vectors and matrices of operators as if the latter were scalars, I m denotes the identity matrix of order m, and i := √ −1 is the imaginary unit. Furthermore, J is a real antisymmetric matrix of order m (we denote the subspace of such matrices by A m ) which specifies the canonical commutation relations (CCRs) for the constituent field processes W 1 , . . . , W m : with J spanning the one-dimensional subspace A 2 of antisymmetric (2 × 2)-matrices, which is an incremental form of the two-point CCRs where the commutator [α, β] := αβ − βα of linear operators α, β is extended to the commutator matrix [ξ, η T ] := ([ξ j , η k ]) 1 j r,1 k s = ξη T − (ηξ T ) T for vectors ξ := (ξ j ) 1 j r , η := (η k ) 1 k s of operators ξ 1 , . . . , ξ r , η 1 , . . . , η s . These CCRs are closely related to the continuous tensor-product structure of the Fock space [44] and are complemented by the commutativity between the Ito increments of W and adapted processes ζ := (ζ(t)) t 0 taken at the same (or an earlier) moment of time: The adaptedness of quantum processes on the system-field space H ⊗ F is understood with respect to a filtration (H t ) t 0 , where and (F t ) t 0 is the Fock space filtration, such that for any t 0 the operators W j (t) act effectively on F t , while X k (t) act on the subspace H t .
The energetics of the quantum system and its interaction with the external fields is specified by a system Hamiltonian H(t) and system-field coupling operators L 1 (t), . . . , L m (t) which are time-varying self-adjoint operators organized as deterministic functions (for ex-ample, polynomials with constant coefficients) of the system variables X 1 (t), . . . , X n (t) and assembled into a vector L(t) := (L k (t)) 1 k m . Accordingly, the operators H(0), L 1 (0), . . . , L m (0) act on the initial space H. Depending on the context, a system operator σ (a function of the initial system variables X 1 (0), . . . , X n (0)) on the initial space H will be identified with its extension σ ⊗ I F to the system-field space H ⊗ F .
The system and the fields form a composite quantum system which evolves according to a quantum stochastic flow as described below. This evolution is specified at any time t 0 by a unitary operator U(t) on H ⊗ F governed by the following stochastic Schrödinger equation [23,24]: with initial condition U 0 := I H⊗F , and so U(t) captures the internal dynamics of the system and the system-field interaction over the time interval In addition, the units are chosen such that the reduced Planck constant is = 1. The quantum stochastic differential equation (QSDE) (6) corresponds to a particular yet important case of the identity scattering matrix in which there is no photon exchange between the fields and the gauge processes [24] can be eliminated from consideration. The term L T 0 dW in (6) can be interpreted as an incremental Hamiltonian of the system-field interaction, while 1 2 L T 0 ΩL 0 dt involves the quantum Ito matrix Ω from (1) and counterbalances the Ito correction term dUdU † = L T 0 ΩL 0 dt (with (·) † denoting the operator adjoint) in the differential relation d( The system variables at time t 0, as operators on the system-field space H ⊗ F , are the images of their initial values under the quantum stochastic flow j t , which maps a system operator σ 0 on the initial space H to the operator on H t in (5). The resulting quantum adapted process σ satisfies the following Hudson-Parthasarathy QSDE [23,24]: where use is made of the Hamiltonian and the coupling operators evolved by the flow j t from (8) as Note that the flow acts on vectors and matrices of operators in an entry-wise fashion. Furthermore, D in (9) is the decoherence superoperator, which acts on σ(t) as The second equality in (11) is applicable to the case where σ is a vector of operators on which the superoperator D acts entry-wise. The superoperator G in (9) is the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) generator [45,46], which is a quantum counterpart of the infinitesimal generators of classical Markov processes [26]. The identity (σ 0 ⊗ I F )U = Uσ, which holds for system operators σ in view of (8) and the unitarity of U(t), allows the QSDE (6) to be represented in the Heisenberg picture by using (10) along with the commutativity (4) as In application to the vector X(t) of system variables, the quantum stochastic flow j t acts entry-wise as in accordance with (7) and (8), and so the corresponding QSDE (9) can be represented in vector-matrix form: where the n-dimensional drift vector F = (G(X j )) 1 j n and the dispersion (n × m) matrix The interaction of the system with the input field W produces the output fields where the system-field unitary evolution is applied to the current input field variables, reflecting the innovation nature of the quantum Wiener process and the continuous tensor product structure of the Fock space mentioned above. The output field Y satisfies the QSDE where the matrix J is given by (2) and L is the vector of the system-field coupling operators from (10). The system-field interaction makes the output field Y different from the input field W only through the drift vector 2JL in (16). The common unitary evolution in (13) and (15) preserves the commutativity between the system variables and the output field variables over the course of time t 0: where the entries of X 0 , W(t) commute as operators on different spaces H, F . For a similar reason, the output field Y inherits the CCR matrix J from the input quantum Wiener process W: where (3) is used with s = t. By means of (1) and (2), the relation (18) can also be established as a corollary of the property that Y in (16) inherits the quantum Ito matrix Ω from W as dYdY T = dWdW T = Ωdt.
In view of (6) (see also (12)), the quantum stochastic flow j t in (8) depends on the system Hamiltonian H 0 and the system-field coupling operators in L 0 . In turn, H 0 and L 0 are usually functions of the initial system variables X 1 (0), . . . , X n (0), such as polynomials or Weyl quantization integrals [27] (see also Equation (32) in [47]). The dependence of H 0 , L 0 on X(0) is inherited by H, L as functions of X and, along with a given commutation structure of the system variables, specifies a particular form of the resulting QSDEs (14) and (16), thereby influencing their tractability.

Open Quantum Harmonic Oscillators with Parametric Dependence
An important class of quantum stochastic systems is provided by multimode open quantum harmonic oscillators (OQHOs) [18], which have an even number n of system variables X 1 , . . . , X n (for example, consisting of n 2 conjugate position-momentum pairs [28]) satisfying the Heisenberg CCRs on a dense domain in H ⊗ F , where the CCR matrix Θ ∈ A n is identified with Θ ⊗ I H⊗F and remains unchanged over the course of time t 0. The Hamiltonian H of the OQHO is a quadratic function and the system-field coupling operators L 1 , . . . , L m in (10) are linear functions of the system variables: which are parameterized by a real symmetric energy matrix R of order n (we denote the subspace of such matrices by S n ) and a coupling matrix N ∈ R m×n . In this case, the QSDEs (14) and (16) take the form where the drifts F = AX and 2JL = CX depend linearly on the system variables X 1 , . . . , X n , with A ∈ R n×n and C ∈ R m×n and the dispersion matrices G = B ∈ R n×m and I m being constant real matrices. This linearity makes several of the dynamic properties of (21) similar to those of a classical linear stochastic system with a state-space realization quadruple (A, B, C, I m ) and the corresponding C m×m -valued transfer function F on the complex plane: allowing for application of transfer function techniques [48]. However, in addition to the noncommutative nature of quantum variables, the matrices of coefficients of these QSDEs have a specific parameterization in terms of the energy, coupling, and CCR matrices: Because the energy matrix R is symmetric and the CCR matrices Θ and J are antisymmetric, the matrices A, B and C satisfy These equalities pertain to the fulfillment of the CCRs (19) and the commutativity (17) at any moment of time, and provide necessary and sufficient conditions for the physical realizability (PR) [21,22] of the linear QSDEs (21) as an OQHO with the CCR matrix Θ for the system variables. The first equality in (24) has the structure of an algebraic Lyapunov equation (ALE) with respect to Θ, which has a unique solution if and only if the Kronecker sum A ⊕ A := I n ⊗ A + A ⊗ I n is a nonsingular matrix, that is, no two eigenvalues of A are centrally symmetric about the origin in C. The latter condition holds, for example, when the matrix A is Hurwitz. For more general open quantum systems (such as anharmonic oscillators, for which the dynamic variables satisfy the CCRs (19) even though the Hamiltonian and the coupling operators are not necessarily quadratic and linear functions of the system variables), the CCR preservation is secured by the unitary evolution of the system variables in (13). By analogy with the state-space realizations of transfer functions in classical linear systems theory [36,37], we use the input-output map for the OQHO (21) with the matrix quadruple (A, B, C, I m ): Now, supposing that the energy and coupling matrices R and N, which specify the Hamiltonian and the coupling operators in (20), depend smoothly on an auxiliary parameter ∈ R (while the quantum Wiener process W is independent of ), then so do the matrices A, B, and C in (23) and the system and output variables which comprise the vectors X(t) and Y(t). Here, the differentiability of X, Y is understood in the weak sense. The corresponding partial derivatives at any t 0 give rise to adapted quantum processes with self-adjoint operator-valued entries satisfying the QSDEs the second of which is in fact an ODE (Y ) = C X + CX involving the time derivative( ) of Y , with zero initial conditions X 0 = 0, Y 0 = 0, because X 0 and Y 0 do not depend on . Here, are the derivatives of the matrices A, B and C from (23) with respect to the parameter satisfying the relations which are obtained by differentiating the PR conditions (24) while the CCR matrices Θ and J remain constant. Because the matrices Θ and J are antisymmetric, whereby By assembling the system variables and their parametric derivatives from (26) to an augmented vector a combination of the QSDEs (21) and (27) allows the parametric derivative of the map (25) to be represented as the input-output map Here, in view of X 0 = 0, the initial condition S 0 = X 0 0 is identified with X 0 , and the matrices A ∈ R 2n×2n , B ∈ R 2n×m and C ∈ R m×2n are given by The same can be obtained by using the transfer functions of the corresponding linear systems (including (22)) as for any s ∈ C which is not an eigenvalue of A. Here, use is made of a particular case of the block matrix inverse formula [49]: The block lower triangular structure of the dynamics matrix A in (32) and (34) is closely related to the Gateaux derivative of the matrix exponential e A in the direction A (see for example [50]): At any time t 0, both X(t) and Y(t) depend linearly on the matrices R and N from (28) in view of the following representation of the solutions of the QSDEs (27), which regards A Xdt + B dW as a forcing term in the first of these QSDEs: where is the unperturbed solution of the first QSDE in (21) which does not depend on R and N . In the case of linear QSDEs with coefficients that depend smoothly on parameters, the mean square differentiability of their solutions with respect to those parameters can be verified directly by using the closed form (37) under certain integrability conditions for the underlying quantum state (specified by a density operator ρ on the system-field space H ⊗ F ) in terms of relevant moments of the system variables such as E(X T 0 X 0 ) < +∞. Here, Eξ := Tr(ρξ) is the quantum expectation which extends entry-wise to vectors and matrices of quantum variables.
An additional insight into the structure of the process X is provided by its crosscommutation relations with X and W. In particular, because X inherits adaptedness from the underlying system variables, it commutes with the Ito increments of the input field W in accordance with (4), and so does S in (31): Furthermore, the differentiation of (19) in , combined with the identity [ξ, η T ] = −[η, ξ T ] T for vectors ξ and η of operators, leads to whereby the cross-commutation matrix [X , X T ] is symmetric, which holds regardless of a particular structure of the QSDEs (21) and remains valid for quantum anharmonic oscillators with nonlinear dynamics. This matrix evolves in time (with zero initial condition, as X 0 = 0) and has a steady-state value computed below. (23) is Hurwitz. Then, there exist the following limits

Theorem 1. Suppose that the matrix A of the OQHO (21) in
where the matrices Θ 21 ∈ S n and Θ 22 ∈ A n are found as unique solutions of the ALEs Proof. In view of the CCRs (19), the commutator matrix for the vector S in (31) is organized as where the blocks Ξ 21 := [X , X T ] = −Ξ T 12 , Ξ 12 := [X, X T ] and Ξ 22 := [X , X T ] consist of time-varying skew self-adjoint operators on the system-field space H ⊗ F . By using the quantum Ito lemma [23,24] and the bilinearity of the commutator along with the QSDE (33) and the commutativity (38), it follows that (43) satisfies the QSDE which, similarly to Equation (56) in [21], reduces to the ODĖ with the initial condition in view of X 0 = 0. Hence, at any time t 0 the matrix Ξ(t) is an imaginary antisymmetric matrix which can be represented as with time-varying matrices Θ 12 (t) = −Θ 21 (t) T ∈ R n×n and Θ 22 (t) ∈ A n . Substitution of (46) into (44) followed by using (34) leads to the Lyapunov ODEṡ with zero initial conditions Θ 21 (0) = 0 and Θ 22 (0) = 0 in accordance with the corresponding blocks in (45). Now, the matrix Θ 21 (t) is symmetric at any time t 0 due to (39). This follows from (47) in view of the symmetry (30). Hence, Θ 12 = −Θ 21 , and the ODE (48) takes the formΘ If the matrix A is Hurwitz, then the solutions of the ODEs (47) and (49) converge to their unique steady-state values which (slightly abusing notation) satisfy the ALEs (41) and (42) and specify the limits (40).
The relations (28)-(36) and Theorem 1 provide infinitesimal perturbation analysis for sensitivity of the internal and output variables of the OQHO to the matrices R and N. Under the perturbations of R and N, the Hamiltonian H and the coupling operators in L given by (20) remain in the corresponding classes of quadratic and linear functions of the system variables. Accordingly, the above analysis is not applicable to more general perturbations, for example, higher order polynomials of the system variables, and is restricted to linear QSDEs, meaning that an alternative approach is needed in the general case.

Transverse Hamiltonian
For the general quantum stochastic system described in Section 2 and governed by (14) and (16), we will now consider its response to arbitrary perturbations in the system Hamiltonian H 0 and the system-field coupling operators in L 0 . More precisely, suppose that they are perturbed in directions K 0 and M 0 as Here, K 0 and the entries of the m-dimensional vector M 0 are self-adjoint operators on the initial system space H, while is a small real-valued parameter as before, and hence, K 0 = H 0 and M 0 = L 0 . In what follows, the derivative (·) := ∂ (·) is taken at = 0. The perturbations K 0 and M 0 in (50) are assumed to be functions of the initial system variables X 1 (0), . . . , X n (0), and these dependencies describe as functions of X(t) under the unperturbed flow (8). For example, in the case of OQHOs considered in Section 3, the perturbations K and M, which are caused by the perturbations in the energy and coupling matrices R and N, inherit the structure of the Hamiltonian as a quadratic function and the coupling operators as linear functions of the system variables in (20), respectively: Returning to the general case, at this stage we avoid technical assumptions regarding K 0 and M 0 in (50); thus, the calculations carried out below for arbitrary perturbations should be regarded as formal. Because the operators H 0 and L 0 completely specify the dynamics of the unitary operator U(t) in (6), which determines the evolution of the system and output field variables, the response of the latter to the perturbations (50) of H 0 and L 0 reduces to that of U(t). Therefore, the propagation of the initial perturbations K 0 and M 0 of the operators H 0 and L 0 through the subsequent unitary system-field evolution can be described in terms of the operator which satisfies V 0 = 0, as U 0 = I H⊗F does not depend on . The smoothness of dependence of U(t) on the parameter is analogous to the corresponding property of solutions of classical SDEs (under suitable regularity conditions [26,51] for their drift and dispersion) and holds at least in the case of linear QSDEs discussed in Section 3. The following theorem is closely related to Stone's theorem on generators of one-parameter unitary groups [52].

Theorem 2.
For any time t 0, the operator V(t) in (53), associated with the unitary evolution U(t) from (6), can be represented as Here, Q(t) is a self-adjoint operator on the system-field space H ⊗ F , which satisfies the zero initial condition Q 0 = 0 and is governed by the QSDE where the imaginary part Im(·) is extended to operators and matrices of operators as Imµ := 1 2i (µ − µ # ), with (·) # denoting the entry-wise operator adjoint. Furthermore, Q(t) depends linearly on the initial perturbations K 0 and M 0 of the Hamiltonian and coupling operators in (50) through their unperturbed evolutions in (51).
Proof. The differentiation of both sides of the unitarity relation U(t) † U(t) = I H⊗F with respect to the parameter at = 0 leads to V † U + U † V = (U † V) † + U † V = 0, which implies self-adjointness of the operator establishing (54). Now, consider the time evolution of Q(t). To this end, the differentiation of (6) with respect to yields where the real part Re(·) is extended to operators and matrices of operators as Reµ := 1 2 (µ + µ # ). By left-multiplying both sides of (57) by U † and recalling (56), it follows that where use is made of the evolved perturbations in (51). By a similar reasoning, a combination of (12) with (56) leads to where use is also made of (58) along with the quantum Ito product rules [24], (1) and (4). It now follows from (56)  In view of Theorem 2, for any fixed but otherwise arbitrary time t 0, the relation (54), represented as has the structure of isolated quantum dynamics in fictitious time , where Q(t) plays the role of a Hamiltonian pertaining to the perturbation of the unitary operator U(t). In order to reflect this property, we will refer to the time-varying operator Q as the transverse Hamiltonian associated with the perturbations K and M of the system Hamiltonian H and the system-field coupling operators in L. The computation of Q is illustrated by the following two examples.

Example 1
In the absence of perturbation of the system-field coupling, when the vector M 0 in (50) consists of zero operators (and hence, so does M in (51)), the transverse Hamiltonian in (61) reduces to Moreover, if the system is isolated, that is, L = 0, then (6) reduces to the ODEU(t) = −iH 0 U(t), which leads to with the Hamiltonian being preserved in time: H(t) = H 0 for all t 0. In this case of isolated system dynamics, (62) takes the form where ad H 0 (·) := [H 0 , ·]. Here, use is made of Hadamard's lemma [28] along with an entire function E(z) := e z −1 z = ∑ +∞ k=0 z k (k+1)! of a complex variable (with E(0) = 1 by continuity) which plays a role in the solution of nonhomogeneous linear ODEs with constant coefficients and constant forcing terms [50]. The Gateaux derivative (53) of (63) can be represented using an operator version of (35) as which provides an alternative verification of (54) for this particular case, with Q given by (64).

Example 2
For the OQHO of Section 3, substitution of the perturbations from (52) into (61) leads to the transverse Hamiltonian which depends linearly on the matrices R and N as well as, in a quadratic fashion, on the past history of the system variables. The latter are given by the unperturbed Equation (37). The last term NΘ, N t= −Tr(ΘN T N )t in (66) (with α, β := Tr(α T β) denoting the Frobenius inner product for real matrices) comes from the relation Im(i(N T ΩN − N T ΩN)) = N T N − N T N, which follows from (1), and the identity X T ΥX = i Υ, Θ , which holds for any matrix Υ ∈ A n in view of the CCRs (19).

Infinitesimal Perturbation Analysis of System Operators
Because the transverse Hamiltonian Q(t) in (56), (61) encodes the propagation of the initial perturbations of the Hamiltonian and coupling operators in (50) through the unitary system-field evolution over the time interval [0, t], it provides a tool for infinitesimal perturbation analysis of general system operators. The following theorem is concerned with an extended setting which, in addition to (50), allows for smooth dependence of a system operator σ 0 on the parameter , meaning that an appropriate infinitesimal perturbation in it is specified by an operator σ 0 on the initial space H, with σ 0 being a function of the initial system variables X 1 (0), . . . , X n (0). Theorem 3. For any self-adjoint system operator σ 0 on the initial space H, which smoothly depends on the same scalar parameter as in (50) and is evolved by the flow (8), the derivative of its evolved version σ(t) with respect to can be represented as at any time t 0, where Q(t) is the transverse Hamiltonian from Theorem 2. Here, the operator φ(t) satisfies the QSDE with zero initial condition φ 0 = 0, where G is the unperturbed GKSL generator from (9) and χ is a linear superoperator given by Proof. By using the Leibniz product rule together with (8), (53), (54) and the self-adjointness of Q(t), it follows that which establishes (67), with the initial condition ϕ 0 = 0 being inherited by φ from Q 0 = 0. We will now obtain a QSDE for the process φ. To this end, by combining the quantum Ito lemma with the bilinearity of the commutator (similarly to the proof of Theorem 1) and using the QSDEs (9) and (55) Here, the quantum Ito product rules are applied together with the commutativity (4) between adapted processes and the Ito increments of the quantum Wiener process W. In the second and third of the equalities (71), use is also made of the relations [α, β T dW] = [α, β T ]dW and [α T dW, β T dW] = 2iIm(α T Ωβ)dt for appropriately dimensioned adapted processes α and β with self-adjoint operator-valued entries. These relations are combined with the identity Im(i[σ, L T ]ΩM) = Re([σ, L T ]ΩM) in the fourth equality of (71). The QSDE (68) is now obtained my multiplying both sides of (71) by i and using the superoperator χ from (69).
As can be seen from (70), the process j t (σ 0 ) in (67) is the response of σ(t) to the initial perturbation σ 0 of the system operator under the unperturbed flow (8) and satisfies the QSDE (9): In contrast to j t (σ 0 ), the operator φ(t) in (67) describes the response of the flow j t itself to the perturbations (50) of the Hamiltonian and coupling operators and will be referred to as the derivative process for the system operator σ. Accordingly, the term i[Q, G(σ)] in the drift of the QSDE (68) is the derivative process for G(σ). At any given time t 0, the last equality in (67) is organized as the right-hand side of the Heisenberg ODE (in fictitious time ) of an isolated quantum system with the state space H ⊗ F and the Hamiltonian Q(t). While j t (σ 0 ) (evolved by the unperturbed flow) depends linearly on σ 0 , the derivative process φ(t) depends linearly on the perturbations K 0 and M 0 of the Hamiltonian H 0 and the coupling operators in L 0 through the transverse Hamiltonian Q and the superoperator χ in (69). In application to the vectors X and L of the system variables and the system-field coupling operators, the QSDEs (68) where F and G are the unperturbed drift vector and the dispersion matrix from (14). In (73), use is also made of the absence X 0 = 0 of initial perturbations in the system variables, whereby j t (X 0 ) = 0 for any t 0. Furthermore, in (74) we have used the relations j t (L 0 ) = j t (M 0 ) = M(t) in view of (50) and (51). Because the matrix J and the quantum Wiener process W do not depend on the parameter , the derivative of the output field Y of the system in (16) evolves according to the ODE where L is governed by the QSDE (74). In particular, for the OQHO from Section 3 with the perturbations in (52), the QSDEs (73)-(75) lead to the relations (27) and (28) which were obtained in Section 3 using more elementary techniques. However, the latter are limited to quadratic perturbations of the Hamiltonian and linear perturbations of the coupling operators in (52), whereas the transverse Hamiltonian approach allows the system response to be investigated for general perturbations of these operators. Therefore, this approach can be used for the development of optimality conditions in quantum control and filtering problems for larger classes of controllers and observers.

Sensitivity of Infinite-Horizon Quantum Averaged Functionals
Similarly to optimal control of classical time invariant stochastic systems [36,37], suppose that the infinite-horizon performance of the quantum system under consideration is described by a cost functional which (whenever it exists) leads to the same Cesaro limit lim t→+∞ ( 1 t t 0 EZ(s)ds) = Z and is to be minimized in optimal quantum control settings. Here, the quantum expectation E(·) is taken over the density operator which is the tensor product of the initial quantum state of the system on the space H and the vacuum state [24] υ on the Fock space F for the external bosonic fields. This expectation is applied in (76) to a quantum criterion process Z, specified as using a function f : R n+m → R. The latter is extended to the noncommutative system variables and the coupling operators (with L being in a bijective correspondence with the drift vector 2JL of the output field in (16), as det J = 0 in view of (2)) to make Z(t) a selfadjoint operator for any t 0. Such an operator-valued extension of f is straightforward in the case of polynomials and can be carried out through the Weyl quantization [27] for more general functions. In the coherent quantum linear-quadratic Gaussian (CQLQG) control and filtering problems [20,35,38], the function f in (78) is a positive semi-definite quadratic form, in which case the minimization of (76) provides an infinite-horizon mean square optimality criterion. Now, if the Hamiltonian and coupling operators of the quantum system are perturbed according to (50), then application of the transverse Hamiltonian Q from Theorems 2 and 3 leads to where φ is the derivative process for Z. Note that despite the equality X 0 = 0, the operator Z 0 can be nonzero due to the dependence of Z 0 = f (X 0 , L 0 ) in (78) on L 0 which is being perturbed. Assuming the existence and interchangeability of appropriate limits, (79) leads to the following perturbation formula for the cost functional Z in (76): The right-hand side of (80) is a linear functional of the perturbations K 0 and M 0 , which describes the corresponding (formal) Gateaux derivative of Z in the direction (K 0 , M 0 ). Therefore, the quantum system is a stationary point of the performance criterion (76) with respect to a subspace T of perturbations (K 0 , M 0 ) in (50) if T is contained by the null space of the linear functional Z in (80): This inclusion provides a first-order necessary condition of optimality in the quantum control problem of minimizing (76) over a manifold of the Hamiltonian and coupling operators with the local tangent space T . While lim t→+∞ Ej t (Z 0 ) in (80) reduces to averaging over the invariant quantum state of the unperturbed system (provided certain integrability conditions are satisfied together with the existence of and weak convergence to the invariant state), the computation of lim t→+∞ Eφ(t) is less straightforward. Due to the product structure (77) of the systemfield state ρ (with the external fields being in the vacuum state υ), the martingale part where ψ is the derivative process for G(Z) and the superoperator χ from (69) is applied to the criterion process Z in (78). The relation (82) is a complicated integro-differential equation (IDE). Nevertheless, this IDE admits an efficient solution, for example, in the case when the system is an OQHO, and the function f in (78) is a polynomial. In this case, due to the structure of the GKSL generator of the OQHO, G(Z) is a polynomial in the system variables of the same degree, leading to a linear relation (with constant coefficients) between the derivative processes φ and ψ and to algebraic closedness in the IDE (82). Therefore, for stable OQHOs (that is, with a Hurwitz matrix A) and a polynomial criterion process Z, the computation of Z and verification of the stationarity condition (81) reduce to averaging over the unperturbed invariant state, which is unique and Gaussian [34]. These calculations are exemplified below.

Example 3
Consider the OQHO from Section 3 described by (19)-(23) with a Hurwitz matrix A. Suppose the criterion process Z in (78) is a quadratic polynomial of the system variables X 1 , . . . , X n and the system-field coupling operators L 1 , . . . , L m : where Π ∈ S + n+m is a given weighting matrix partitioned into blocks Π 11 ∈ S + n , Π 22 ∈ S + m and Π 12 = Π T 21 ∈ R n×m , and is an auxiliary matrix which involves the coupling matrix N from (20). Here, S + r denotes the set of real positive semi-definite symmetric matrices of order r. Then, in view of the relation X 0 = 0, it follows from (83) and (52) that The second equality in (83) allows the quantum average of the corresponding derivative process φ in (79) to be represented as Here, Υ is the expectation of the derivative process i[Q, Ξ] for Ξ and, as such, satisfies the following IDE, similar to (82):Υ where the superoperator χ from (69) is applied to Ξ entry-wise as with Θ • denoting the th row of the CCR matrix Θ. By substituting the Hamiltonian and coupling operators of the OQHO from (20) into the GKSL generator G in (9) and (11) and using the CCRs (19) together with the state-space matrices (23), it follows that Substitution of (89) into (87) reduces the IDE to a Lyapunov ODE: Because the matrix A is Hurwitz, the unperturbed OQHO has an invariant state which is Gaussian [34] with zero mean EX = 0 and a complex positive semi-definite Hermitian covariance matrix E(XX T ) = Σ + iΘ of order n, for which the real part Σ ∈ S + n is a unique solution of the ALE AΣ + ΣA T + BB T = 0.
Therefore, under appropriate integrability conditions for χ(Ξ), the solution of the Lyapunov ODE (90) has a limit Υ ∞ := lim t→+∞ Υ(t), which is a unique solution of the ALE where lim t→+∞ Eχ(Ξ) can be computed by averaging χ(Ξ) over the invariant Gaussian state of the OQHO. Therefore, by assembling (85) and (86) into (80), it follows that where the limit reduces to averaging over the invariant Gaussian state and is expressed in terms of E(XM T ) as The relations (84), (88) and (92)-(94) allow the Gateaux derivative Z of the quadratic cost functional (76), specified by (83), to be computed through mixed moments of the system variables X and the perturbations K and M over the invariant zero-mean Gaussian quantum state, the covariance matrix of which can be found from (91). In particular, if K and M are polynomials of the system variables X 1 , . . . , X n , then these moments can be calculated in terms of the covariances by using the Isserlis-Wick theorem [53]. Alternatively, the perturbations K and M can be trigonometric polynomials, that is, linear combinations of unitary Weyl operators [27] associated with the system variables, or more generally represented as the Weyl quantization integrals similar to those in Equation (32) of [47]. Here, α and β are countably additive measures of finite variation on the σ-algebra B n of Borel subsets of R n , which take values in C and C m , respectively, and satisfy the Hermitian property α(S) = α(−S) and β(S) = β(−S) (where (·) is the complex conjugate) for any S ∈ B n , ensuring that K and the entries of M in (96) are self-adjoint operators in view of the second equality in (95). Now, the Weyl CCRs W u+v = e iu T Θv W u W v for all u, v ∈ R n (with their infinitesimal Heisenberg form given by (19)) imply that ∂ u W u = ∂ v e iv T (X+Θu) v=0 W u = i(X + Θu)W u , whereby Because the quantum expectation commutes with the differential operator on the righthand side of (97) and the quasi-characteristic function [54] of the invariant zero-mean Gaussian state is given by which can also be obtained using quantum Price's theorem [55]. A combination of the second equality from (96) with (98) leads to The computation of the limit in (93) in the Weyl quantization framework can now be completed by substituting (99) into (94). In this framework, the matrix Eχ(Ξ), which is needed for finding Υ ∞ from the ALE (92) for the last term in (93), is computed by averaging (88) in a similar fashion.

Mean Square Optimal Coherent Quantum Filtering Problem
We will now demonstrate an application of the transverse Hamiltonian variational approach from Sections 4-6 to a quantum filtering problem for the open quantum system described in Section 2, referred to as a quantum plant. The plant has the Hamiltonian H, vector L of system-field coupling operators, input field W, vector X of internal variables, and output field Y governed by the QSDEs (14) and (16). Although the plant is not necessarily an OQHO, we assume that the plant variables X 1 , . . . , X n are organized as position-momentum pairs and satisfy the CCRs (19).
Suppose that the quantum plant is cascaded in a measurement-free field-mediated fashion with another open quantum system, which plays the role of a coherent quantum observer and is driven by the plant output Y and another quantum Wiener process ω of an even dimension µ on a different symmetric Fock space F; see Figure 1. The series connection of a quantum plant with a coherent quantum observer, mediated by the plant output field Y and affected by the environment through the input quantum Wiener processes W and ω. The observer design objective is that the drift part of its output η has to approximate the plant variables of interest in a mean square optimal fashion.
The observer is endowed with its own initial Hilbert space H, and dynamic variables ξ 1 (t), . . . , ξ ν (t) with a CCR matrix Λ ∈ A ν such that [ξ(t), ξ(t) T ] = 2iΛ for any t 0, similarly to (19), as well as an (m + µ)-dimensional output field η(t): As the observer is driven by the plant output Y together with the quantum noise ω, over the course of time the observer output η acquires quantum statistical correlations with the plant variables and can be used for estimating the latter in a mean square optimal fashion as specified below. To this end, we denote the observer Hamiltonian by Γ and the vectors of operators of coupling of the observer with the plant output Y and the quantum Wiener process ω by respectively. The Hamiltonian Γ and the coupling operators Φ 1 , . . . , Φ m , Ψ 1 , . . . , Ψ µ are functions of the dynamic variables ξ 1 , . . . , ξ ν of the observer and hence (by the commutativity [X, ξ T ] = 0) commute with the plant variables and functions thereof, including H and L.
The plant and observer form a composite quantum stochastic system with a vector X of dynamic variables satisfying the CCRs and driven by an augmented quantum Wiener process W with the following Ito table: dWdW T = Ωdt, Ω : Here, is the Ito matrix of the quantum Wiener process ω of the observer, which is defined similarly to Ω from (1) and (2) as The Ito matrix Ω of the process W in (103) is block-diagonal due to the commutativity [W, ω T ] = 0 and the statistical independence between the quantum Wiener processes W and ω, which act on the different Fock spaces F and F and are assumed to be in the appropriate vacuum states. The quantum feedback network formalism [7] allows the Hamiltonian H of the plant-observer system and its vector L of operators of coupling with W to be computed as Hence, the internal and output variables ξ 1 , . . . , ξ ν and η 1 , . . . , η m+µ of the observer in (100) are governed by the QSDEs Here, is the GKSL generator of the plant-observer system and is the corresponding decoherence superoperator in accordance with (11), (103) and (105). In (106), use is made of the partial decoherence superoperator ∆, which acts on the observer variables as in view of (102)-(105). Now, consider a coherent quantum filtering (CQF) problem formulated as the minimization of the mean square discrepancy between r linear combinations of the plant variables of interest and the entries of the drift part 2J L of the observer output η in (107), as specified by given weighting matrices S ∈ R r×n and T := T 1 T 2 with T 1 ∈ R r×m and T 2 ∈ R r×µ . The criterion process Z in (108) is similar to that in (83): with X being the only subvector of X from (102) which is present in (109). The minimization in (108) is over the observer Hamiltonian Γ and the vector Φ of the observer-plant coupling operators in (101) as functions of the observer variables ξ 1 , . . . , ξ ν . This mean square optimal CQF problem extends [20] in that we do not restrict attention to linear observers even if the plant is an OQHO and Ψ depends linearly on the observer variables ξ. Note that the Hamiltonian H of the plant and its coupling L to the input quantum noise W are fixed, as is the coupling Ψ of the observer to the input quantum noise ω; see Figure 1. If Γ 0 and Φ 0 are perturbed in the directions consisting of self-adjoint operators representable as functions of the initial observer variables ξ 1 (0), . . . , ξ ν (0), then the corresponding perturbations of the plant-observer Hamiltonian and coupling operators in (105) are By applying Theorem 2 and using (103) In particular, if both the plant and the unperturbed observer are OQHOs with Hurwitz dynamics matrices (while the perturbations in (110) are not necessarily linear-quadratic), then Z can be found in a form similar to (93) along the lines of Example 3 from Section 6, which is due to the criterion process Z in (108) and (109) being a quadratic function of the plant and observer variables in X from (102). The transverse Hamiltonian method outlined above was used in [41] to show that in the mean square optimal CQF problem for linear quantum plants those observers which are locally optimal in the class of linear quantum observers cannot be improved locally in the sense of the first-order optimality conditions (81) by varying the Hamiltonian and coupling operators of the observer along linear combinations of the Weyl operators associated with the observer variables.

Concluding Remarks
For a class of open quantum systems governed by Markovian Hudson-Parthasarathy QSDEs, we have introduced a transverse Hamiltonian and derivative processes associated with perturbations of the system Hamiltonian and system-field coupling operators along with the influence of such perturbations on the system dynamics. This provides a fully quantum tool for infinitesimal perturbation analysis of system operators and cost functionals with infinite-horizon formulations that involve averaging over the invariant quantum state of the system. The proposed variational method has been accompanied by detailed proofs and can be used for the development of first-order necessary conditions of optimality in quantum control and filtering problems. We have illustrated these ideas for OQHOs with quadratic performance criteria and the mean square optimal CQF problem. The results of the paper are also applicable to perturbation analysis and variational problems in optimal control and filtering synthesis for more complicated networks of quantum stochastic systems, such as in [7,56].

Acknowledgments:
The author thanks the anonymous reviewers for useful comments.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this paper: