{'fraction_non_alphanumeric': 0.05066004929268375, 'fraction_numerical': 0.024713866099463182, 'mean_word_length': 4.347928139387921, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Invariant states of linear quantum stochastic systems under Weyl perturbations of the Hamiltonian and coupling operators

This paper is concerned with the sensitivity of invariant states in linear quantum stochastic systems with respect to nonlinear perturbations. The system variables are governed by a Markovian Hudson-Parthasarathy quantum stochastic differential equation (QSDE) driven by quantum Wiener processes of external bosonic fields in the vacuum state. The quadratic system Hamiltonian and the linear system-field coupling operators, corresponding to a nominal open quantum harmonic oscillator, are subject to perturbations represented in a Weyl quantization form. Assuming that the nominal linear QSDE has a Hurwitz dynamics matrix and using the Wigner-Moyal phase-space framework, we carry out an infinitesimal perturbation analysis of the quasi-characteristic function for the invariant quantum state of the nonlinear perturbed system. The resulting correction of the invariant states in the spatial frequency domain may find applications to their approximate computation, analysis of relaxation dynamics and non-Gaussian state generation in nonlinear quantum stochastic systems.


Introduction
At the core of emerging quantum computation and quantum information technologies is the employment of hidden resources of physical systems at the atomic scales governed by the laws of quantum mechanics. In contrast to the classical Newtonian mechanics, such systems are described by operator-valued quantities acting on an underlying Hilbert space which evolve according to unitary similarity transformations when the system is isolated from the environment. This unitary evolution is specified by the Hamiltonian which represents self-energy of the system. When the system is coupled to the surroundings in the form of interaction with a macroscopic measuring apparatus, the internal state of the system is disturbed in a random fashion, depending on what is measured, which makes noncommuting operators inaccessible to simultaneous measurement [18,34,51]. This specific stochastic nature of quantum systems is described in terms of quantum probability which replaces the scalar-valued classical probability measures with density operators (quantum states) acting on the same Hilbert space as the dynamic variables.
An important part of quantum computation protocols [39] is the initialization of quantum systems in certain classes of states. Similarly to the classical terminal state control problem [49], the quantum state preparation is a quantum control problem which can be solved by varying the parameters of the Hamiltonian either in an openloop fashion or by using feedback [5]. Another approach [70,71] to the state preparation consists in letting the quantum plant of interest interact with external bosonic fields (such as nonclassical electromagnetic radiation) in a specifically arranged fashion, resulting in a dissipative system with a required invariant state. The interaction of This work is supported by the Air Force Office of Scientific Research (AFOSR) under agreement number FA2386-16-1-4065. arXiv:1711.03503v1 [math-ph] 9 Nov 2017 an open quantum system with such fields is modelled by using the Hudson-Parthasarathy quantum stochastic calculus [22,43]. This approach represents the open quantum dynamics in the form of quantum stochastic differential equations (QSDEs) driven by noncommutative counterparts of the classical Wiener process [27], with the quantum Wiener processes acting on the symmetric Fock spaces [45]. In addition to the Hamiltonian, the energetics of the system-field interaction is specified by the coupling operators. In combination with scattering matrices (which take into account photon exchange between the fields in multichannel setups), this approach is used by the theory of quantum feedback networks [13,24] in order to model interconnections of open quantum systems which interact with each other and the environment.
The quantum state generation problem is important both for finite-level systems (such as qubit registers) and continuous variables systems which find application in different (for example, quantum-optical) platforms of quantum computing [39]. An important class of such systems is provided by open quantum harmonic oscillators (OQHOs) whose dynamic variables satisfy canonical commutation relations (CCRs) similar to those of the quantum mechanical positions and momenta [34,51], and are governed by linear QSDEs. The Hamiltonian and coupling operators of the OQHO are quadratic and linear functions of the system variables. The dynamics of such systems, when they are driven by bosonic fields in the vacuum state, are quite similar (at least at the level of the mean values and covariances) to those of classical Gaussian Markov diffusion processes (for example, the Ornstein-Uhlenbeck process [27]). In particular, if the OQHO is initialized in a Gaussian quantum state [44], then its Gaussian nature is preserved in time. Furthermore, if the dynamics matrix of the linear QSDE is Hurwitz, the state of the system converges weakly [3,4] to a unique invariant Gaussian state regardless of whether the initial state is Gaussian. This property of OQHOs enables them to be used for generating a Gaussian state by allowing the system to evolve over a sufficiently long period of time. The dissipation mechanism, which is built in this procedure, secures stability of the invariant state being achieved in the long run [42].
The implementation of OQHOs (for example, using quantum optical components such as cavities, beam splitters and phase shifters) can involve parametric uncertainties, leading to Gaussian invariant states with different covariances [67]. More complicated inaccuracies result in deviations from the idealized linear-quadratic energetics. The appearance of nonquadratic terms in the Hamiltonian and nonlinearities in the coupling operators makes the governing QSDEs nonlinear and leads to non-Gaussian invariant quantum states. The deviation from the nominal invariant Gaussian quantum state can be studied in the framework of a particular description of the uncertainties in the energy operators. Such uncertainties can be modeled as functions of the system variables. However, there are many ways to extend usual functions of several real or complex variables to the noncommutative quantum variables. One of such extensions is provided by the Weyl functional calculus [8] which employs unitary Weyl operators whose role in this context is similar to that of the spatial Fourier harmonics. In fact, the Weyl quantization replaces these harmonics in the Fourier integral representations of classical functions with the Weyl operators. In application to the energy operators of the quantum stochastic system, the Weyl quantization leads to partial differential equations (PDEs) or integro-differential equations (IDEs) for the evolution of the quasi-characteristic function (QCF) [4]. The latter is a quantum counterpart of the usual characteristic functions of classical probability distributions and encode information on the mixed moments of quantum system variables.
The QCF and its Fourier transform -the Wigner quasi-probability density function (QPDF) -constitute the basis of the Wigner-Moyal phase-space approach [16,38] which avoids the "burden of the Hilbert space" by representing the quantum state dynamics in the more conventional form of PDEs and IDEs involving only real or complex variables and functions thereof. Although the Moyal equations [38] for the QCF and QPDF dynamics were obtained originally for isolated quantum systems, there also are extensions of the phase-space approach to different classes of open quantum systems [14,15,29,33,62].
In the present paper, following [62], we apply the phase-space approach to the infinitesimal perturbation analysis of the invariant state of an OQHO when its linear-quadratic energy operators are subject to variations represented in the Weyl quantization form. The effect of such perturbations of the Hamiltonian on the second moments of the system variables (with the coupling operators remaining unperturbed) was studied in [53]. In the case of linear system-field coupling, the Weyl quantization of the Hamiltonian was used for approximate computation of invariant states via operator splitting [32,54] in [62] and dissipation relations for a weighted L 2 -norm of the QCF, similar to the norm in the Bessel potential space [55], and for the χ 2 -divergence (the second-order Renyi relative entropy [50]) of the actual QPDF from the nominal Gaussian PDF. Similarly to and in furtherance of the previous developments, we make advantage of the tractability of the Gaussian state dynamics for linear quantum stochastic systems and the ease of computing the moments of the Weyl operators over Gaussian quantum states. This is combined with the machinery of linear PDEs and related semigroups [7,56] together with the generalized functions [57] in the spatial frequency domain in order to develop a perturbative (linear response) theory for QCFs of quantum systems subject to the Weyl variations both in the Hamiltonian and coupling operators. More precisely, we compute the influence functions which relate the first-order corrections of the QCF and the first two moments for the perturbed invariant quantum state to the strength functions of the Weyl variations and discuss norm bounds for these perturbation operators.
We also mention that the Weyl quantization of the energy operators was applied to the measurement-based quantum stochastic filtering in [63,64] and was used in [60,61] in regard to optimality conditions for coherent (measurement-free) quantum control and filtering problems for linear quantum plants. The coherent quantum control and filtering by interconnection [25,35,37,41,46,48,47,58,59,66,73] provide a promising modern paradigm which has the potential to outperform the traditional observation-actuation approaches of the classical control theory.
The paper is organised as follows. Section 2 specifies the class of quantum stochastic systems being considered. Section 3 describes the Weyl quantization of energy operators and the integro-differential equation for the evolution of the quasi-characteristic function. Section 4 discusses the dynamics of quantum states and preservation of their Gaussian nature for open quantum harmonic oscillators. Section 5 specifies the Weyl perturbations of the Hamiltonian and coupling operators of such an oscillator. Section 6 computes the corresponding first-order correction terms of the QCF and QPDF for the resulting perturbed system. Section 7 discusses upper bounds for the norms of the linear operators relating these corrections to the variations in the energy operators. Section 8 gives numerical examples of the invariant state perturbation analysis, including the sensitivity bounds. Section 9 provides a conclusion and outlines possible directions of research.

Quantum stochastic dynamics
We will be concerned with an open quantum system interacting with external fields as shown in Fig. 1 dynamics are described in terms of an even number n of variables X 1 , . . . , X n assembled into a vector X := (X k ) 1 k n (vectors are organized as columns). These system variables are time-varying self-adjoint operators on a complex separable Hilbert space H . Their evolution is governed by a Markovian Hudson-Parthasarathy quantum stochastic differential equation (QSDE) [22,43] whose structure is described below (we will often omit the time arguments for brevity). Although it formally resembles classical SDEs [27], the QSDE (1) is driven by a vector W := (W k ) 1 k m of an even number m of self-adjoint operator-valued quantum Wiener processes W 1 , . . . ,W m acting on a symmetric Fock space F . These processes model the external bosonic fields [17,43] and satisfy the quantum Ito relations where ⊗ is the Kronecker product of matrices, and I r is the identity matrix of order r. In application to matrices of operators, the transpose (·) T acts as if their entries were scalars. In contrast to the identity diffusion matrix of the standard Wiener process, Ω is a complex positive semi-definite Hermitian matrix with an orthogonal antisymmetric imaginary part ImΩ = J (so that J 2 = −I m ). Therefore, the quantum Wiener processes W 1 , . . . ,W m do not commute with each other. Moreover, they have the following two-point commutator matrix where [α, β ] := αβ − β α is the commutator of linear operators α and β . In accordance with W representing the m-channel input field, h := (h k ) 1 k m in (1) (1) consists of n operators and is obtained by the entrywise application of the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) generator [12,30], which acts on a single system operator ξ as and involves the system Hamiltonian h 0 (which is also a self-adjoint operator on H ) in addition to the coupling operators h 1 , . . . , h m . The energy operators h 0 , h 1 , . . . , h m are appropriately extended functions of the system variables X 1 , . . . , X n . The superoperator L in (1) is a quantum analogue of the infinitesimal generators of classical Markov diffusion processes [27,55]. The specific structure of the drift vector and the dispersion matrix of the QSDE (1) comes from the system-field interaction which drives a unitary operator U(t) on the system-field tensor-product space H := H 0 ⊗ F (with H 0 the initial space for the action of the system variables at time t = 0) as where U(0) = I H is the identity operator on H . The unitary operator U(t), which is associated with the systemfield interaction over the time interval from 0 to t, is adapted in the sense that it acts effectively on the subspace H 0 ⊗ F t , where {F t : t 0} is the Fock space filtration. The corresponding quantum stochastic flow at time t involves a unitary similarity transformation, which acts on operators ζ on H as ζ → U(t) † ζU(t) (and applies entrywise to vectors of operators) and evolves the system variables as where (·) † is the operator adjoint. The QSDE (1) can be obtained from (8) by using (7) and the quantum Ito formula [22,43] in combination with (3) and commutativity between the forward Ito increments dW (t) and adapted processes taken at time s t. Adapted processes ξ , which are functions of the system variables, satisfy QSDEs of the same form An alternative representation of the quantum Wiener process W (which reflects photon absorption and emission and is particularly relevant to quantum optics) is in terms of the field annihilation a 1 , . . . , a m/2 and creation a † 1 , . . . , a † m/2 processes [22,43], assembled into vectors a := (a k ) 1 k m/2 and a # := (a † k ) 1 k m/2 : Here, (·) # denotes the entrywise operator adjoint (reducing to the complex conjugate in the case of complex matrices), and the real and imaginary parts are extended to matrices N with operator-valued entries as ReN = 1 2 (N + N # ) and ImN = 1 2i (N − N # ) which consist of self-adjoint operators. Also, is an auxiliary doubled-up vector. In accordance with (2) and (3), the quantum Ito table of the annihilation and creation processes is represented in terms of (11) as where (·) † := ((·) # ) T applies to matrices of operators as the transpose of the entrywise adjoint and reduces to the complex conjugate transpose (·) * := ((·)) T for complex matrices. Also, the term which is part of the diffusion term in (7), is related to a different vector L := (L k ) 1 k m/2 of (not necessarily self-adjoint) coupling operators L 1 , . . . , L m/2 , so that Here, use is also made of the property J 2 = −I m of the matrix J in (3) mentioned above. The relations (10) and (12) allow the two alternative representations of the external fields and the system-field coupling operators to be used interchangeably. As a result of the joint system-field evolution described by the unitary operator U(t) from (7), the output field Y := (Y k ) 1 k m is given by and satisfies the QSDE dY = 2Jhdt + dW.
Here, the vectors b := (b k ) 1 k m/2 and b # := (b † k ) 1 k m/2 consist of the corresponding output annihilation and creation operators: In view of (10), (12) and (14), the processes b and b # satisfy the QSDEs which are related to each other by conjugation. The QSDEs (1) and (14) describe a particular yet important scenario of quantum stochastic dynamics with the identity scattering matrix [22,43]. Endowed with additional features (for example, more general scattering matrices describing photon exchange between the fields), such QSDEs are employed in a unified formalism [13,24] for modelling feedback networks of quantum systems which interact with each other and the external fields. The unitary evolution in (8) and (13) preserves the commutativity between the system and output field variables in the sense that In view of (5) and (14), the process Y inherits from W the two-point commutator matrix: Despite the system-output commutativity (15), the noncommutativity of the output fields in (17) makes them inaccessible to simultaneous measurement. However, they can be fed as an input to other open quantum systems without conversion of the quantum operator-valued processes to classical real-valued signals. The resulting coherent field-mediated interconnection can be used in fully quantum communication channels in the form of cascades of (possibly distant) quantum systems. Such cascade connections are also employed for the generation of certain classes of Gaussian quantum states (see, for example, [31,67]).

Evolution of the quasi-characteristic function
Whereas the QSDE (9) pertains to the Heisenberg picture of quantum dynamics, its dual Schrödinger picture version, known as the master equation [69], describes the evolution of the reduced density operator. The latter is a quantum counterpart of the probability distribution for a classical Markov diffusion process, which evolves according to the Fokker-Planck-Kolmogorov equation (FPKE) specified by the infinitesimal generator of the process [27,55]. However, in contrast to the classical case, the master equation carries the "burden of the Hilbert space" in the sense that it involves operator-valued quantities rather than usual functions of real variables (such as a probability density function). The Wigner-Moyal phase-space approach [16,38] reconciles these two pictures by dealing with the reduced density operator indirectly through the moments of the system variables. These moments (and their time evolution) can be recovered from the quasi-characteristic function [4,21] considered below. To this end, the system variables are assumed to satisfy the Weyl canonical commutation relations (CCRs) for any time t 0. These are represented in terms of the unitary Weyl operators [8] which inherit their time dependence from the system variables. Here, Θ is a constant nonsingular real antisymmetric matrix of order n which specifies the commutator matrix as an infinitesimal form of the Weyl CCRs (18). Note that (18) is closely related to the Baker-Campbell-Hausdorff formula for the exponentials of more general operators which commute with their commutator (see, for example, [10, pp. 128-129]).
Example 1. Let the system variables consist of pairs of conjugate position and momentum operators [34] comprising the vectors q := (q k ) 1 k n/2 and p := −i∂ q = (−i∂ q k ) 1 k n/2 , with the reduced Planck constant being set to = 1, and q 1 , . . . , q n/2 denoting the Cartesian coordinates in the position space R n/2 . Then, without loss of generality, they can be assembled into the vector (where, as before, the time arguments are omitted for brevity). The corresponding CCR matrix Θ in (18) and (20) takes the form and, up to a factor of 1 2 , coincides with the symplectic structure matrix in classical Hamiltonian systems [2]. In what follows, Eζ = Tr(ρζ ) denotes the quantum expectation over the system-field tensor-product density operator where ϖ is the initial quantum state of the system, and υ is the vacuum state [43] of the input fields. The averaging of the Weyl operator W u , associated with the system variables by (19), leads to the complex-valued quasicharacteristic function (QCF) Here, the second equality describes the Hermitian property of Φ(t, u) with respect to its spatial argument u and follows from the second equality in (19). Also, the relation W 0 = I H implies that Φ(t, 0) = 1. However, in view of the Weyl CCRs (18), the Bochner-Khinchin positiveness criterion [11] for the characteristic functions of classical probability distributions is replaced with its quantum mechanical weighted version [4,21]: the complex Hermitian matrix e iu T j Θ u k Φ(t, u j − u k ) 1 j,k is positive semi-definite for arbitrary points u 1 , . . . , u ∈ R n and any = 1, 2, 3, . . .. The spatial Fourier transform F of (24) yields a real-valued quasi-probability density function (QPDF) Although the function is not necessarily nonnegative everywhere [20] (since, as mentioned above, the QCF Φ itself does not have to be positive semi-definite), it satisfies the normalization condition R n (t, x)dx = Φ(t, 0) = 1 and is a quantum analogue of classical PDFs. In particular, Φ and encode the mixed moments E(X k 1 × . . . × X k ) of the system variables for any 1 k 1 , . . . , k n, provided Φ(t, u) is times continuously differentiable with respect to u ∈ R n . For example, the first and second-order moments of the system variables are expressed in terms of the QCF Φ and the QPDF similarly (modulo taking the real part) to the corresponding moments of classical random variables: where ∂ u (·) and ∂ 2 u (·) denote the gradient vector and the Hessian matrix with respect to u ∈ R n . Furthermore, the QCF Φ in (24) can be used for evaluating the generalized moments which involve nonlinear (but not necessarily polynomial) functions of the system variables represented in the Weyl quantization form. Here, σ : R n → C is a given function, which specifies such a moment and can be a generalized function [57]. For example, (26) and (27) can be obtained from (28) and (29) by letting σ be distributional partial derivatives (of up to the second order) of the n-dimensional Dirac delta function δ .
Since the QCF Φ in (24) is the expectation of the Weyl operators, its time evolution satisfies the integrodifferential equation (IDE) This relation is obtained by averaging the QSDE (9), which is applied to ξ := W u , with the martingale part −i[ξ , h T ]dW not contributing to the right-hand side of (30) (since the forward Ito increments of W commute with adapted processes, as mentioned before, and the input fields are in the vacuum state).
In general, the right-hand side of (30) is not easily related to the QCF Φ. However, (30) becomes an algebraically closed equation for the QCF in the framework of the Weyl quantization model for the energy operators of the system. More precisely, following [62], we assume that the system Hamiltonian h 0 and the system-field coupling operators h 1 , . . . , h m in (6) are obtained by the Weyl quantization [8] of real-valued functions on R n with the Fourier transforms H k : R n → C as where W u is the Weyl operator (19), and arg(·) denotes the argument of a complex number. The second equality in (31) employs the Hermitian property of the functions H k (that is, H k (−u) = H k (u) for all u ∈ R n ), which, in combination with the second equality in (19), ensures that the operators h k are self-adjoint. Accordingly, the vector h of the system-field coupling operators h 1 , . . . , h m is related to the vector-valued map H : If the function H k is absolutely integrable, then (31) can be understood as a Bochner integral [72] which yields a bounded operator h k whose operator norm admits an upper bound due to unitarity of the Weyl operators W u for all u ∈ R n . Polynomial functions of the system variables can be obtained by letting the Fourier transforms H k in (31) be the derivatives of the Dirac delta function (as mentioned before in regard to (26)- (29)).
The following theorem 1 uses an integral operator A which maps a function ϕ : R n → C to the function A(ϕ) : R n → C given by Here, the kernel function V : R n × R n → C is computed as where H 0 and H are the Fourier transforms from (31) and (32), and the function K : R n × R n → R m×m is expressed as in terms of the CCR matrix Θ in (18) (1) have the Weyl quantization form (31). Then the IDE (30) takes the form where the integral operator A, described by (33)-(35), acts over the spatial argument of the QCF Φ.
Proof. By substituting the Weyl quantization of the system Hamiltonian h 0 from (31) into the GKSL generator L in (6), it follows that Here, use is made of the Lie algebraic property of the Weyl operators, which follows from the CCRs (18). In combination with (38), the Weyl quantization of the system-field coupling operators in (32) leads to where (18) is used again. Here, the kernel function V 1 : R n × R n → C is organised as a weighted convolution of the function H with itself: By a similar reasoning, Here, in view of the antisymmetry of the CCR matrix Θ and the Hermitian property of the quantum Ito matrix Ω , the kernel function V 2 : R n × R n → C is computed as From (40) and (42), it follows that where use is made of the function K from (35). A combination of (37), (39) and (41) allows the superoperator L in (6) to be evaluated at the Weyl operator as where is given by (34) in view of (43). The averaging of (44) represents (30) (36), with the integral operator A in (33).
If the system and fields are uncoupled, the IDE (36) reduces to the Moyal equation for isolated quantum systems [38]: which follows from (33) by letting H = 0 in (34). Note that (36) and its special case (45) are phase-space representations of the master equation for the reduced density operator in the spatial frequency domain. The corresponding IDE for the QPDF in (25) is obtained from (36) by a unitary similarity transformation of the operator A (due to unitarity of the Fourier transform F up to a factor of (2π) n/2 ): where F −1 is the inverse Fourier transform over the spatial variables in R n . The kernel function Π of the integral operator F is related to V in (33) by

Linear-Gaussian dynamics of open quantum harmonic oscillators
We will now consider the case when the Hamiltonian h 0 is a quadratic function and the coupling operators h 1 , . . . , h m are linear functions of the system variables described by where R is a real symmetric matrix of order n, and M ∈ R m×n . The matrices R and M will be referred to as the energy and coupling matrices, respectively. The corresponding Fourier transforms H 0 and H in (31) and (32) are given by where δ and δ are the distributional gradient vector and Hessian matrix of the n-dimensional Dirac delta function δ , and N, L := Tr(N * L) denotes the Frobenius inner product [19] of real or complex matrices.
Example 2. Suppose the system variables are the quantum mechanical positions and momenta as in (21) of Example 1, and the system Hamiltonian consists of the quadratic potential and kinetic energy parts as with a stiffness matrix K and a positive definite mass matrix M. Then the representation (48) holds with a blockdiagonal energy matrix Regardless of a particular structure of the system variables and the energy matrix, considered in the above example, the combination of the CCRs (20) and the linear-quadratic energetics in (48) and (49) makes the QSDEs (1) and (14) linear: where A ∈ R n×n , B ∈ R n×m and C ∈ R m×n are constant matrices of coefficients: The special dependence of A, B, C on the energy and coupling matrices R and M from (48) and (49), the symmetry of R and the antisymmetry of the matrices Θ and J, lead to the physical realizability (PR) conditions [25,52]: the first of which is closely related to the preservation in time of the CCRs (20), while (60) pertains to (15). Such quantum stochastic systems are referred to as open quantum harmonic oscillators (OQHOs). Their QCF and QPDF dynamics can be obtained from Theorem 1 by using the Weyl quantization and the generalized functions in the spatial frequency domain as the following lemma shows for completeness.
Lemma 1 Suppose the initial system variables of the OQHO, described by (54)- (58) and driven by vacuum input fields, have finite second moments (that is, E(X(0) T X(0)) < +∞). Then the QCF Φ(t, ·) in (24) remains twice continuously differentiable with respect to its spatial variables for any time t 0 and satisfies the first-order linear PDE Also, the QPDF in (25) is governed by the parabolic PDE where div(·) is the divergence operator over the spatial variables (which acts on matrix-valued functions f := ( f jk ) 1 j m,1 k n on R n with m > 1 in a row-wise fashion as div f := ∑ n k=1 ∇ k f jk 1 j m , where ∇ k is the partial derivative with respect to the kth Cartesian coordinate in R n ).
Proof. Substitution of (50) and (51) into (34) allows the action of the integral operator A in (33) on a twice continuously differentiable function ϕ to be represented in the form where div v (·) and div s (·) are the divergence operators with respect to the vectors v, s ∈ R n , and use is made of the matrices A and B from (56) and (57) together with the symmetry of R and antisymmetry of Θ . Here, we have also used the relation for a continuously differentiable function f := ( f jk ) 1 j m,1 k n : R n → C m×n and the distributional gradient of the n-dimensional Dirac delta function (see, for example, [57]). Now, if the initial system variables of the OQHO have finite second moments, then the linear QSDE (54) propagates this mean square integrability to all subsequent times as where · is the operator norm of a matrix, and the factor m = TrΩ comes from the quantum Ito matrix in (3). This follows from the relation the product structure (23) of the system-field density operator and the input fields being in the vacuum state. The property E(X(t) T X(t)) < +∞ in (65) secures the twice continuous differentiability of the QCF Φ(t, ·) with respect to its spatial argument for any t 0. Therefore, in view of (63), the IDE (36) takes the form of the PDE (61). Accordingly, the IDE (46) for the QPDF with the kernel function (47) can be reduced to the PDE (62) by applying the spatial Fourier transform to both sides of (61).
Note that (62) coincides with the FPKE for a classical Markov diffusion process ξ in R n governed by the SDE dξ = Aξ dt + Bdω (where ω is the standard Wiener process in R m ) with the infinitesimal generator F † which acts on a smooth function ϕ : R n → R (with the gradient ϕ and the Hessian matrix ϕ ) as (see, for example, [27,55]). The PDEs (61) and (62) show that the OQHO inherits the preservation of the Gaussian nature of quantum states from classical Markov processes governed by linear SDEs. The QCF of a Gaussian quantum state [44] is parameterized by the mean vector µ ∈ R n and the real part Σ of the quantum covariance matrix as where u Σ := | √ Σ u|. The Gaussian QCF (68) is identical to its classical counterpart except that Σ satisfies the generalized form Σ + iΘ 0 of the Heisenberg uncertainty principle [18], which is a stronger property than the positive semi-definiteness of Σ . The corresponding Gaussian QPDF is given by provided Σ 0 (the latter condition also makes Φ µ,Σ square integrable). Therefore, if the system is initialized at a Gaussian quantum state, then the state remains Gaussian over the course of time, with the parameters of the QCF (68) satisfying the ODEsμ where the matrices A and B are given by (56) and (57). Furthermore, in this case, the system variables and the output fields, described by (13)- (17), are in an augmented Gaussian state at every moment of time. Furthermore, if A is Hurwitz, then the system has a unique invariant quantum state. This state is Gaussian with the mean vector µ = 0 (in view of (70)) and the real part of the quantum covariance matrix being the unique steady-state solution of the Lyapunov ODE in (71): which is the infinite-horizon Gramian of the matrix pair (A, B), satisfying the algebraic Lyapunov equation (ALE) The property that the matrix P is physically meaningful (that is, satisfies the Heisenberg uncertainty relation) can be verified directly by noting that in view of (59), (72) and the positive semi-definiteness of the quantum Ito matrix Ω . The convergence to the invariant Gaussian quantum state also holds in the sense of weak convergence of probability measures [3] (appropriately modified for the quantum mechanical case [4]) irrespective of whether the initial state is Gaussian. The invariant Gaussian QCF and QPDF (with (76) assuming that (A, B) is controllable which ensures that P 0) are the fixed points of the corresponding semi-groups associated with the IDEs (36) and (46) in the sense that for all t 0. In the case of OQHOs being considered, the operator exponential e tA in (77) can be computed by solving the first-order PDE (61) with an initial condition ϕ : R n → C through the method of characteristics [7,56] as Here, use is made of (68) together with the finite-horizon controllability Gramian of the pair (A, B) satisfying the ODE (71) with zero initial condition Σ (0) = 0: with TrΣ (t) being bounded by the second term on the right-hand side of (65). The Gramians Σ (t) and P = lim t→+∞ Σ (t) in (72) satisfy the identity e tA Pe tA T + Σ (t) = P, t 0.
The operator exponential e tA in (78) can also be obtained directly from (66) (regardless of whether the system variables have finite second moments) as Here, use is made of the commutativity between the input fields (more precisely, their increments) and the initial system variables, the product structure of the system-field state (23) and the quasi-characteristic functional of the vacuum fields: for any t 0 and any locally square integrable function f : R + → R m . Accordingly, e tF is an integral operator whose kernel function is the fundamental solution [7,56] of the FPKE (62), provided (A, B) is controllable, which is equivalent to Σ (t) 0 for all t > 0. This operator maps an initial condition ψ : R n → R of the FPKE to the function e tF (ψ)(x) = R n e tA y,Σ (t) (x)ψ(y)dy, where use is made of the Gaussian PDF (69). The function y → e tA y,Σ (t) (x) is the Markov transition kernel of the classical diffusion process with the infinitesimal generator (67). Therefore, any matrix pair (A, B) satisfying (72) and the PR condition (59), with A Hurwitz, can be used in order to generate a zero mean Gaussian state with a given quantum covariance matrix P + iΘ in (74) as an invariant state of the OQHO governed by (54). Note that inaccuracies in the energy and coupling matrices R and M of such a system will only perturb the matrix P for this invariant state without distorting its Gaussian nature. Their effects are considered, for example, in [67]. However, the Hamiltonian and coupling operators can, in general, be subject to a broader class of perturbations which bring the quantum system beyond the linear-quadratic energetics (48) and (49).

Weyl variations of the energy operators
Now, consider a class of quantum stochastic systems, which are similar to those in [53,62,64] and result from perturbing the linear-quadratic energy operators (48) and (49) of the nominal OQHO of the previous section subject to the following variations: Here, ε is a small real-valued parameter, and h 0 and h := ( h j ) 1 j m describe nonquadratic and nonlinear parts of the perturbed energy operators h ε 0 , h ε 1 , . . . , h ε m , which are assumed to be represented in the Weyl quantization form (see Section 3) as These Weyl variations depend on d n system variables comprising a subvector SX (where S ∈ {0, 1} d×n is a suarray of a permutation matrix of order n) and are determined by the Fourier transforms Ψ : R d → C and ϒ := (ϒ k ) 1 k m : R d → C m of functions on R d with values in R and R m , respectively. In addition to being Hermitian, both functions Ψ and ϒ are assumed to be absolutely integrable, so that (84) and (85) are Bochner integrals which produce self-adjoint bounded operators on the system-field Hilbert space H . The corresponding Fourier transforms in (31) and (32) are given by where Since SS T = I d , the matrix S T describes an isometry between R d and the subspace S T R d ⊂ R n . Accordingly, the integrals in (88) and (89), as generalized functions [57], are measures on this subspace with densities Ψ and ϒ (with respect to the d-dimensional Lebesgue measure on S T R d ) and values in C and C m , respectively. Regardless of the particular form (84) and (85) of the perturbations h 0 and h in (82) and (83), the perturbed GKSL generator (6) depends on ε in a quadratic fashion: where L := L 0 is the unperturbed generator. Here, the parametric derivative L := ∂ ε L ε ε=0 is a linear superoperator which acts on a system operator ξ as where M is the coupling matrix of the nominal OQHO in (49). In the case of the Weyl variations (84) and (85), the superoperator L depends linearly on the functions Ψ and ϒ which quantify the strength of such perturbations of the energy operators.
Example 3. Suppose ϒ = 0 (that is, the coupling operators in (83) remain unperturbed). Then the QSDE (1) takes the form [62] where the nonlinear dependence of the drift term on the system variables comes from the nonquadratic part of the Hamiltonian in (82). In the particular case (21) and (52), considered in Examples 1 and 2, the Weyl variation of the potential energy leads to the Hamiltonian which corresponds to (82) with the energy matrix R given by (53) and the matrix S extracting the position variables from the vector X: This perturbation enters the QSDE (92) through the momentum part in the form of an additional "force" term in the drift vector as where Here, use is made of the CCR matrix Θ from (22), whereby 2Θ S T = 0 −I d , and the matrices A and B in (56) and (57) are partitioned into blocks A jk ∈ R d×d and B j ∈ R d×m as Note that the form of the position QSDE (96) remains unperturbed. The integrals in (93) and (98) lend themselves to closed-form evalution, for example, if Ψ is a linear combination of quadratic-exponential functions which describe Gaussian-shaped "wells" or "bumps" in the potential energy [9] (see also [62,Section 9]). More precisely, suppose the function Ψ is given by where α is a real-valued constant, γ ∈ R d , and Λ is a real positive definite symmetric matrix of order (95). The corresponding integral in (93) is the inverse Fourier transform of (100): whose right-hand side (up to a constant factor) is the Gaussian PDF with the mean vector γ and covariance matrix Λ . The parameter γ of the potential φ in (101) specifies the location of a centre of attraction (if α < 0) or repulsion (if α > 0) in the position space R d with the stiffness matrix φ (γ) = −αΛ −1 . The exponentially fast decay of this potential (together with its gradient) at infinity resembles the Morse potential [36]. The corresponding perturbation force term in (97) is obtained by a straightforward differentiation of (101) with respect to q as φ (q) = φ (q)Λ −1 (γ − q). The nonquadratic perturbations of the Hamiltonian, considered in the above example, are applicable to the modelling of open quantum systems with multiextremum potential energy landscapes in the form which extends (101). In addition to being ubiquitous in macromolecular dynamics simulations, multiextremum potentials are also of interest from the viewpoint of adiabatic quantum computing [1], where they are represented by the terminal Hamiltonian of a quantum system, towards which the initial Hamiltonian is adiabatically evolved together with its ground state. In fact, the efficiency of the quantum annealing protocols in comparison with the classical optimization algorithms is particularly important for the minimization of complicated functions of many variables.

First-order correction of the quantum state
In accordance with the perturbed energy operators in (82) and (83), the operators A, F in (33), (46), the QCF Φ and the QPDF acquire dependence on the parameter ε and will be denoted by A ε , F ε , Φ ε , ε , respectively (the subscript (·) 0 for the unperturbed quantities at ε = 0 will be omitted for brevity). Furthermore, in view of (34) and (35), the integral operators A ε , F ε depend on ε in a quadratic fashion, inheriting such dependence from the GKSL generator L ε in (90). The QCF and QPDF of the perturbed quantum system admit the asymptotic series expansions where the functions satisfy the recurrent IDEs for all r = 0, 1, 2, . . .. These are obtained by applying the Leibniz product rule to the differentiation of the equations (36) and (46) with respect to ε, starting from the PDEs (61) and (62) for the QCF Φ (0) = Φ and the QPDF (0) = of the unperturbed oscillator. Here, ( r ) k := r! k!(r−k)! is the binomial coefficient, and use is made of the parametric derivatives which vanish for all r > 2 in view of the quadratic dependence of the operators A ε and F ε on ε, thus leading to the band structure of the equations (103) and (104). This structure can be represented symbolically in vector-matrix form as (3) (4) . . .
In what follows, we will consider only the first-order infinitesimal perturbations of the quantum state in (102) described by Their evolution is governed by with zero initial conditions Φ(0, ·) = 0 and (0, ·) = 0, where Φ and are the QCF and QPDF for the unperturbed OQHO satisfying the PDEs (61) and (62). Here, the operators A and F from (105) act on the functions (106), and their parametric derivatives are related by the same unitary similarity transformation as in (46). The operator B is related to the superoperator L in (91) through the QCF Φ in (24) and the Weyl operator W u in (19) by The following lemma establishes the structure of B as a perturbation of the integral operator A in (33) rather than through the averaging in (110).
Lemma 2 Suppose the Hamiltonian and coupling operators of the OQHO are subject to the Weyl variations (82)-(85). Then B in (109) is an integro-differential operator which acts on a continuously differentiable function ϕ : R n → C as Proof. The kernel function of the operator B can be obtained by substituting (86) and (87) into (34) and differentiating the result with respect to ε: where the function K is given by (35) and H 0 and H are the measures in (88) and (89). The integrals in (113) originate from differentiating the integral in (112). The integral operator B with the kernel function V acts on a function ϕ : R n → C as Here, the integral operators B 1 and B 2 are associated with the Weyl variations of the system-field coupling operators as where the computations are similar to (63) and also employ the relation (64) together with the structure of the function K in (35). From (115) and (116), it follows that In combination with the structure of the functions H 0 and H as the measures in (88) and (89), substitution of (117) into (114) leads to (111) .
Application of the Duhamel formula [7,55] to the solutions of the initial value problems for the linear evolutionary equations (107) and (108) (with zero initial conditions) leads to for any time t 0. The linear operators on the right-hand sides of (118) and (119) also involve the parametric derivatives of the exponential maps [28] from (78) to the Gaussian QCF and QPDF Φ * and * of the nominal invariant state described in Section 4, with the limits being understood similarly to the weak convergence [3], so that the convergence in (120) is pointwise. Here, −A −1 B and −F −1 G denote the operators which relate the steady-state solutions of (107) and (108) to Φ * and * subject to the normalization constraints Since the QPDF perturbation * in (121) is related by the spatial Fourier transform F to the QCF perturbation Φ * , we will be concerned mainly with the latter. This will also allow advantage to be made of the fact that the operator exponential e tA , associated with the fundamental solution of the first-order linear PDE (61), is simpler than e tF for the second-order PDE (62).
Theorem 2 Suppose the unperturbed OQHO has a Hurwitz matrix A in (56), and its Hamiltonian and coupling operators are subject to the Weyl variations (82)-(85). Then the corresponding perturbation (120) of the invariant QCF is related by to the strength functions Ψ and ϒ of the Weyl variations. Here, the functions F : R n ×R d → R and G := (G k ) 1 k m : where Here, the function K is given by (35), M is the coupling matrix for the unperturbed oscillator in (49), and P is the infinite-horizon controllability Gramian from (72).
Proof. By combining the second equality in (120) with the structure of the operator exponential e tA in (78), it follows that where use is made of the finite-horizon controllability Gramian Σ (t) from (79). Evaluation of the operator B from (111) at the unperturbed invariant Gaussian QCF Φ * in (75) leads to where use is also made of the relation Φ * (v) = −Φ * (v)Pv for all v ∈ R n (following from the quadratic-exponential structure of Φ * ) and the matrix D from (126). Substitution of (128) into (127) establishes the integral representation (123) with the kernel functions By using (75) and (80), it follows that for any t 0 and u, v ∈ R n . Substitution of (131) with v := S T w into (129) and (130), leads to (124) and (125), where the convergence of the integrals is secured by the exponentially fast decay of the integrands since the matrix A is Hurwitz.
The functions F and G in (124) and (125) depend on the parameters of the unperturbed oscillator and can be identified with the Frechet derivatives ∂ Ψ Φ * and ∂ ϒ Φ * of the QCF perturbation Φ * with respect to the strength functions Ψ and ϒ of the Weyl variations in the energy operators. These derivatives describe the response of the invariant state of the quantum system in the framework of the infinitesimal perturbation analysis and allow the QCF corrections to be calculated for particular strength functions. Example 4. In the setting of Example 3 (with unperturbed system-field coupling operators), we will now compute the QCF correction in response to the Gaussian-shaped variation (101) of the potential energy (see also Examples 1 and 2). By substituting its Fourier transform Ψ from (100) together with the influence function F from (124) into (123) and letting ϒ = 0, it follows that Here, use is made of the moment-generating function for a zero-mean Gaussian distribution in R d with the covariance matrix Ξ −1 , where is a real positive definite symmetric matrix of order d. In view of the structure of the matrix S in (94), the block P 11 ∈ R d×d of the matrix P := P 11 P 12 from (72) is associated with the position variables in accordance with the partitioning (21) and (99). Also, σ + , σ − : R + × R n → C d are auxiliary functions given by σ ± (t, u) := S(P ± iΘ )e tA T u + iγ where (134) is combined with the structure of the CCR matrix Θ from (22). Since the functions in (135) converge to a common limit lim t→+∞ σ ± (t, u) = σ ± (0, 0) = iγ exponentially fast due to A being Hurwitz, the integral in (132) is convergent for any u ∈ R n . The presence of double exponentials (such as e −u T e tA PS T w ) in the integrals (124), (125) and (132) complicates their closed-form evaluation. Upper bounds for the influence functions F and G, which relate the QCF perturbation to the Weyl variations in Theorem 2, will be discussed in Section 7. These functions also encode more tractable first-order corrections to the mean vector and the second-moment matrix for the system variables in the invariant (in general, non-Gaussian) state under such variations: where E ε (·) denotes the quantum expectation over the invariant state of the perturbed system. Note that the constant imaginary part ImE ε (XX T ) = Θ does not contribute to the parametric derivative in (136). Therefore, the corresponding truncations of the asymptotic series expansions for these moments in the perturbed invariant state take the form E ε X = ε µ + . . . , ReE ε (XX T ) = P + ε P + . . . , where use is made of the zero mean vector and the quantum covariance matrix P + iΘ of the unperturbed invariant Gaussian state. In view of the first of the relations (137), the system variables in the perturbed invariant state acquire nonzero mean values.
Here, the matrix-valued functions Q : R d → R n×n and R 1 , . . . , R m : R d → R n×n are, respectively, quadratic and cubic polynomials defined as the unique solutions of the ALEs for all w ∈ R d , where S(N) := 1 2 (N + N T ) denotes the symmetrizer of square matrices, the matrix D is given by (126), and (·) k• is the kth row of a matrix.
Proof. By differentiating the representation (123) for the QCF perturbation Φ * in accordance with (26), (27) and (136), it follows that where ∂ u G(u, w) is the Jacobian (m × n)-matrix of G, and ∂ 2 u G k (u, w) is the Hessian matrix of the kth coordinate function of G with respect to u ∈ R n . Now, (124) leads to where (similarly to (72) and (73)) the last integral is the solution Q(w) of the ALE (140) which depends parametrically on w ∈ R d in a quadratic fashion. Also, (125) yields for all k = 1, . . . , m, where A −T := (A −1 ) T , and the matrix R k (w) is the solution of the ALE (141) which is a cubic function of w ∈ R d . Substitution of (144) and (146) into (142) leads to (138). Similarly, substitution of (145) and (147) into (143) establishes (139).
Note that in the setting of Example 3 (where only the Hamiltonian is perturbed), the mean vector correction µ can be obtained by averaging the drift of the QSDE (92) over the unperturbed invariant Gaussian state with the QCF Φ * in (75) as in accordance with (138) with ϒ = 0. This is yet another demonstration of the origin of the quadratic-exponential factors in the correction terms.

Norm bounds for the perturbation operators
For what follows, we will now specify a class of the strength functions Ψ and ϒ of the Weyl variations (84) and (85) in the energy operators. More precisely, they are assumed to be square integrable with a weight in the sense that |||Ψ ||| θ < +∞, |||ϒ ||| θ < +∞ for some θ > 0. The parameter θ specifies the following weighted L 2 -norm [62] for a scalar or vector-valued function ϕ on R d : The fulfillment of (148) ensures that Ψ and ϒ are not only square integrable, but are also absolutely integrable in view of the relations with , which follow from (149) and the Cauchy-Bunyakovsky-Schwarz inequality. Note that the limiting case ||| · ||| 0 is the standard L 2 -norm · 2 . Now, the relations (142) and (143) in the proof of Theorem 3 show that the functions (144)-(147) specify the Frechet derivatives of the first-order corrections µ and P for the mean vector and the covariance matrix with respect to Ψ and ϒ . In particular, the function w → i∂ u F(u, w) u=0 , computed in (144), describes the Frechet derivative ∂ Ψ µ with respect to square integrable strength functions Ψ of the Weyl variations (84) in the system Hamiltonian. Similarly to (150), the sensitivity of the mean vector of the perturbed invariant quantum state to Ψ with the norm |||Ψ ||| θ can be quantified by Here, use is made of an auxiliary zero-mean Gaussian PDF f with the covariance matrix 1 2 (θ I d + SPS T ) −1 given by Since the matrix S is of full row rank, the controllability of the matrix pair (A, B) in (56) and (57) (whereby P in (72) satisfies P 0) ensures that SPS T 0. In this case, the right-hand side of (151) has a limit as θ → 0, which corresponds to the standard L 2 -norm ∂ Ψ µ 2 . In view of (145)-(147), the norms of the other Frechet derivatives ∂ Ψ P, ∂ ϒ µ, ∂ ϒ P are related in a similar fashion to the higher-order moments for the Gaussian PDF (152), which can be computed in closed form with the aid of the Isserlis theorem [23,26]. The details of these calculations are cumbersome and are omitted for brevity (see also [65]). Beyond the first two moments of the perturbed invariant quantum state discussed above, its sensitivity to the Weyl variations satisfying (148) can be quantified for the QCF correction Φ * by the corresponding L 2 -induced norm of the integral operator (Ψ ,ϒ ) → Φ * in (123) of Theorem 2. Since this norm is hard to compute, we will consider an upper bound which uses a weighted Hilbert-Schmidt norm [6] of this operator involving the weighted L 2 -norms of the influence functions F and G in (124) and (125): Here, with a slight abuse of notation, the weighting is over the second arguments of the functions F and G. Up to a factor of (2π) −n/2 , the right-hand side of (153) is also an upper bound for the norm * 2 = (2π) −n/2 Φ * 2 of the first-order correction (121) of the invariant QPDF (due to the Plancherel identity applied to * as the Fourier transform of Φ * ). Although the evaluation of the integrals (124) and (125), required for (154) and (155), can be carried out numerically, we will consider upper bounds which give a qualitative insight into the influence functions F and G. To this end (similarly to [68,Theorem 6]), let λ and Γ be a positive scalar and a real positive definite symmetric matrix of order n satisfying the algebraic Lyapunov inequality Such pairs (λ ,Γ ) exist since the matrix A is Hurwitz. Moreover, (156) has a positive definite solution Γ for any λ < − max 1 k n Reλ k , where λ 1 , . . . , λ n are the eigenvalues of A. In a generic case when A is diagonalizable, a particular choice of such a matrix Γ can be carried out by using the eigenbasis of A.
Lemma 3 Suppose the unperturbed OQHO has a controllable matrix pair (A, B) in (56) and (57), with A Hurwitz, and its Hamiltonian and coupling operators are subject to the Weyl variations (82)-(85). Then the influence functions F and G in (124) and (125) admit the upper bounds for all u ∈ R n and w ∈ R d . Here, (λ ,Γ ) is any pair of a positive scalar λ and a real positive definite symmetric matrix Γ of order n satisfying (156), and where r(·) denotes the spectral radius of a square matrix.
Proof. For any pair (λ ,Γ ), satisfying (156), application of the Gronwall-Bellman lemma argument in the form of the differential matrix inequality e 2λt e tA Γ e tA T = e 2λt e tA (AΓ +Γ A T + 2λΓ )e tA T 0 leads to e tA Γ e tA T e −2λt Γ for all t 0, which is equivalent to the contraction in terms of the matrix operator norm. By combining (124) with submultiplicativity of the operator norm and using (160), it follows that for all u ∈ R n and w ∈ R d , thus establishing (157), where the quantity is given by (159). Here, the positive definiteness of the matrix SPΓ −1 PS T is ensured by S being of full row rank and by P 0 (due to controllability of the pair (A, B)) together with Γ 0. Note that the double exponential function in the second last integral of (161) comes from (160) and is circumvented by changing the integration variable. By a similar reasoning, a combination of (125) with (160) leads to the upper bound which proves (158). In (162), we have used the property J = 1 of the matrix J from (4), whereby the function K in (35) satisfies K(u, v) 2 for all u, v ∈ R n .
Since the upper bounds in (157) and (158) involve the quadratic-exponential factors which are proportional to Gaussian PDFs, they lead to upper bounds for the weighted L 2 -norms of the influence functions F and G in (154) and (155) which can be evaluated numerically by Monte-Carlo methods. More precisely, a combination of (154) with (157) leads to Here, the expectation is over independent zero-mean Gaussian random vectors ν and ω with values in R n and R d and covariance matrices 1 2 P −1 and 1 2 (θ I d + SPS T ) −1 , respectively (so that ω has the PDF f given by (152)). By a similar reasoning, (155) and (158) imply that |||G||| 2 −θ 2 λ 2 π (n+d)/2 det P det(θ I d + SPS T ) Therefore, the right-hand sides of (163) and (164) can be estimated via the Monte-Carlo method by sampling from the classical Gaussian distributions.

NUMERICAL EXAMPLES OF INVARIANT STATE PERTURBATION ANALYSIS
As an illustration of the above results, we will now provide two numerical examples of infinitesimal perturbation analysis of invariant quantum states under the Weyl variations of the energy operators. First, consider a one-mode OQHO with n = 2 system variables consisting of the position q and momentum p operators (see Example 1) driven by m = 2 external fields. The energy and coupling matrices R and M of the OQHO in (48) and (49) were randomly generated so as to make the matrix A in (56) The stiffness of this well (at its centre γ) is −αΛ −1 = 919.0101. The perturbed potential 1 2 Kq 2 + εφ (q) (with ε = 0.03) is shown in Fig. 2. The real and imaginary parts of the corresponding QCF correction Φ * , evaluated by numerical integration of (132) using (133)-(135) with Ξ = 2.3797, are depicted in Fig. 3. The QPDF correction * , obtained by the Fourier transform of Φ * in (121), is shown in Fig. 4. The corresponding marginal QPDF correction +∞ −∞ * (x 1 , x 2 )dx 2 for the position variable of the perturbed OQHO is depicted in Fig. 5. Although the perturbed marginal QPDF of the position variable is no longer Gaussian, it remains nonnegative everywhere on the real line (as does the classical probability distribution of any self-adjoint operator). The correction in Fig. 5 increases the value of the perturbed QPDF in the vicinity of the secondary potential well in Fig. 2 at x 1 ≈ 3 by "taking" some probability from negative values of the position coordinate x 1 (in order to satisfy the normalization constraint similar to (122)). Also, this correction contributes a positive shift +∞ −∞ x 1 * (x 1 , x 2 )dx 1 dx 2 = 1.3677 to the mean value of the position from 0 towards the well.    As another example, consider a two-mode OQHO with n = 4 system variables consisting of the conjugate positions q 1 , q 2 and momenta p 1 , p 2 driven by m = 4 fields. The energy and coupling matrices of the OQHO (also randomly generated so as to make the matrix A in (56)  is {−1.7804 ± 5.1767i, −1.2680 ± 0.8111i}. Now, let the potential energy part of the system Hamiltonian be subject to the Weyl variations whose strength functions Ψ satisfy (148). With the matrix S given by (94) with d = 2, the norm |||∂ Ψ µ||| −θ in (151), which quantifies the sensitivity of the perturbed invariant mean vector, is shown in Fig. 6.

Conclusion
For a class of open quantum harmonic oscillators, whose linear-quadratic coupling and Hamiltonian operators are subject to the Weyl variations leading to nonlinearities in the governing QSDE and non-Gaussian state dynamics, we have obtained the first-order correction terms for the QCF and QPDF of the invariant state. We have also discussed the norms of the linear operators which relate these corrections to the Weyl variations. These infinitesimal perturbation analysis results can find applications to the Gaussian state generation and approximate invariant state computation taking into account inaccuracies and uncertainties in the implementation of such systems. Another possible direction of the research, reported in this paper, is an extension of these ideas to more general states of the external fields (such as nonvacuum Gaussian states [40]).