Applied Mathematics and Optimization manuscript No.
(will be inserted by the editor)

Multi-point Gaussian States, Quadratic-exponential Cost Functionals, and Large
Deviations Estimates for Linear Quantum Stochastic Systems∗

arXiv:1707.09302v1 [math.OC] 28 Jul 2017

Igor G. Vladimirov · Ian R. Petersen · Matthew R. James

Abstract This paper is concerned with risk-sensitive performance analysis for linear quantum stochastic systems interacting with external bosonic fields. We consider a cost functional in the form of the exponential moment of the
integral of a quadratic polynomial of the system variables over a bounded time interval. An integro-differential equation is obtained for the time evolution of this quadratic-exponential functional, which is compared with the original
quantum risk-sensitive performance criterion employed previously for measurement-based quantum control and filtering problems. Using multi-point Gaussian quantum states for the past history of the system variables and their first
four moments, we discuss a quartic approximation of the cost functional and its infinite-horizon asymptotic behaviour.
The computation of the asymptotic growth rate of this approximation is reduced to solving two algebraic Lyapunov
equations. We also outline further approximations of the cost functional, based on higher-order cumulants and their
growth rates, together with large deviations estimates. For comparison, an auxiliary classical Gaussian Markov diffusion process is considered in a complex Euclidean space which reproduces the quantum system variables at the level of
covariances but has different higher-order moments relevant to the risk-sensitive criteria. The results of the paper are
also demonstrated by a numerical example and may find applications to coherent quantum risk-sensitive control problems, where the plant and controller form a fully quantum closed-loop system, and other settings with nonquadratic
cost functionals.
Keywords Linear quantum stochastic system · Gaussian quantum state · Risk-sensitive quantum control
Mathematics Subject Classification (2010) 81S25 · 81S05 · 81S22 · 81P16 · 81P40 · 81Q93 · 81Q10 · 60G15

1 Introduction
The main theme of the present paper is a class of risk-sensitive performance criteria for linear quantum stochastic
systems. Such systems, also referred to as open quantum harmonic oscillators (OQHOs) [15], play the role of building
blocks in linear quantum systems theory [53]. This paradigm is part of the broader area of quantum filtering and control
(see, for example, [3, 4, 7, 8, 13, 27, 28, 30, 46, 66]) which is concerned with achieving certain dynamic properties for
open quantum systems interacting with surroundings such as classical measuring devices, other quantum systems or
external quantum fields. In particular, these properties may include stability, optimality and robustness. In contrast
to their classical counterparts, quantum systems are equipped with noncommuting operator-valued variables whose
evolution obeys the laws of quantum mechanics [35, 42, 55] and statistical characteristics are described in terms of
quantum probability [43]. The applications include, for example, artificially engineered systems for influencing the
state of matter at atomic scales through its interaction with nonclassical light in quantum optics [63] and quantum
computing [45].
∗ This work is supported by the Air Force Office of Scientific Research (AFOSR) under agreement number FA2386-16-1-4065. The authors are with

the Research School of Engineering, College of Engineering and Computer Science, Australian National University, Canberra, Acton, ACT 2601,
Australia.
Igor G. Vladimirov
E-mail: igor.g.vladimirov@gmail.com
Ian R. Petersen
E-mail: i.r.petersen@gmail.com
Matthew R. James
E-mail: matthew.james@anu.edu.au

2

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

A unified language for the modelling of open quantum systems, interacting with the environment, is provided
by Hudson-Parthasarathy quantum stochastic differential equations (QSDEs) [23, 47, 51] (see, also, the review paper
[19] and references therein) which govern the system variables in the Heisenberg picture of quantum dynamics. The
QSDEs are driven by quantum Wiener processes on a symmetric Fock space [47, 49] which represent external bosonic
fields (such as quantised electromagnetic radiation). The structure of these quantum stochastic dynamics depends on
the energetics of the system-field interaction and the self-energy of the system. These are captured by the systemfield coupling operators and the system Hamiltonian, with the latter specifying the internal dynamics of the system
in isolation from the environment. This approach to open quantum systems and their interconnections is widely used
in quantum control [3, 13, 28, 30], where it is also combined with the theory of quantum feedback networks [18, 29]
(see also [69] and references therein). In the case of OQHOs, the system variables satisfy canonical commutation
relations (similar to those for the quantum-mechanical position and momentum operators [42]), the Hamiltonian is a
quadratic function and the coupling operators are linear functions of the variables. The linearity of the resulting QSDEs
makes them tractable in many respects including the dynamics of the first two moments of the system variables and
preservation of the Gaussian nature of system states (provided the fields are in the vacuum state). Due to the linearGaussian quantum dynamics, such systems resemble Gaussian Markov diffusion processes generated by classical
linear SDEs [32]. This analogy is exploited in the quantum counterparts [13, 38, 46] of the classical linear quadratic
Gaussian (LQG) control and filtering problems [2, 34] (and also H∞ -control settings [30]). The optimality in these
approaches is understood as the minimization of mean square cost functionals which are concerned with second-order
moments of the system variables at one instant in time (or the integrals of such moments in finite-horizon formulations
[60]). A more general class of cost functionals for quantum systems, which use one-point averaging of nonquadratic
functions of system variables, is considered in [52].
Qualitatively different are the quantum risk-sensitive performance criteria [27, 28] (see also [11, 67]) which were
used previously for measurement-based quantum control and filtering problems. They employ weighted mean square
values of time-ordered exponentials satisfying operator differential equations (where noncommutativity of quantum
variables plays its part). Similarly to their classical predecessors [5, 26, 64] and in contrast to the mean square values
of the system variables themselves, the quantum risk-sensitive cost functionals involve multi-point quantum states and
higher-order moments for the system variables at different instants. In addition to being a challenge from a purely
theoretical point of view, the computation and minimization of such functionals is also of practical importantance.
In fact, this approach leads to controllers and filters which secure a more conservative behaviour of the system not
only in terms of the one-point second-order moments of the system variables but also their higher-order multi-point
moments. Also, there are connections (though more limited in the quantum case) between the risk-sensitive criteria
and robustness with respect to statistical uncertainty in the driving noise described in terms of relative entropy (see,
for example, [27, 67] and references therein).
A quadratic-exponential functional (QEF), considered in the present paper, is organised as the exponential moment
of the integral of a quadratic form of the system variables over a bounded time interval. It differs from the original
quantum risk-sensitive criterion [27, 28], mentioned above, in that a weighted mean square value of the time-ordered
exponential is replaced with the mean value of the operator exponential of the integral. At the same time, the cost
functional studied here is, in fact, a more straightforward extension of its classical predecessors [5, 26, 64] to the
quantum case and also imposes penalty both on the second and higher-order moments of the system variables. As
before, the shift towards the higher-order moments is controlled by a risk-sensitivity parameter. These features motivate
the study of such cost functionals for quantum systems.
To this end, we obtain an integro-differential equation for the time evolution of the QEF and compare it with the
original quantum risk-sensitive performance criterion. Assuming that the OQHO is driven by vacuum fields and using
the multi-point Gaussian quantum states for the system variables at different instants and their first four moments, we
study a quartic approximation of the cost functional and its infinite-horizon asymptotic behaviour. The computation
of the asymptotic growth rate for this approximation is reduced to solving two algebraic Lyapunov equations. This
is a quantum extension of similar results of [59] on classical linear stochastic systems. However, we also discuss an
auxiliary classical Gaussian Markov diffusion process in a complex Euclidean space which reproduces the quantum
system variables at the covariance function level but has different higher-order moments relevant to the risk-sensitive
criteria. In addition to the quartic approximation of the QEF, we also outline its higher-order approximations and their
asymptotic growth rates along with a large deviations estimate of the Cramer type [12, 57].
The paper is organised as follows. Section 2 specifies the class of linear quantum stochastic systems being considered and provides background material. Section 3 discusses multi-point Gaussian quantum states associated with the
system variables at different instants. Section 4 defines the quadratic-exponential functional and establishes an integrodifferential equation for its time evolution. Section 5 discusses a quartic approximation of the QEF and its state-space
computation. Section 6 outlines further approximations of the QEF using higher-order cumulants. Section 7 establishes
an upper bound for the cumulants and provides a related large deviations estimate. Section 8 considers a correspon-

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

3

dence between vectors of self-adjoint quantum variables and classical random vectors at the level of covariances and
discusses its violation in regard to higher-order moments. Section 9 provides a numerical example which demonstrates
the quartic approximation of the QEF for a two-mode oscillator. Section 10 makes concluding remarks. Appendices A,
B and C provide auxiliary lemmas on computing the commutator and covariance of quadratic functions of Gaussian
quantum variables, and also on averaging in a class of convolution-like integrals of matrix-valued functions.

2 Linear quantum stochastic systems
We consider a quantum system whose dynamic variables X1 , . . . , Xn are time-varying self-adjoint operators, defined
on a dense domain in a complex separable Hilbert space H and satisfying the Weyl canonical commutation relations
(CCRs)
T
Wu Wv = eiv Θ u Wu+v
(1)
√
for all u, v ∈ Rn . Here, i := −1 is the imaginary unit, and Θ := (θ jk )16 j,k6n is a nonsingular real antisymmetric
matrix of order n (which is assumed to be even). Also,
Wu := eiu

TX

(2)

is the unitary Weyl operator [14] parameterised by a vector u := (uk )16k6n ∈ Rn whose entries specify the linear
combination uT X = ∑nk=1 uk Xk of the system variables assembled into the vector
 
X1
 .. 
(3)
X :=  . 
Xn
(vectors are organised as columns unless indicated otherwise). Due to self-adjointness of uT X, the adjoint of the Weyl
operator in (2) satisfies Wu† = W−u . Also, (1) implies that
[Wu , Wv ] = −2i sin(uTΘ v)Wu+v ,

(4)

where [ξ , η] := ξ η − ηξ is the commutator of linear operators. In view of (4), the Heisenberg infinitesimal form of
the Weyl CCRs (1) is described by the commutator matrix
[X, X T ] := ([X j , Xk ])16 j,k6n
= XX T − (XX T )T = 2iΘ ⊗ IH

(5)

(on a dense domain in H), where ⊗ is the tensor product, and IH is the identity operator on H (the subscript will
sometimes be omitted for the sake of brevity). The transpose (·)T acts on matrices with operator-valued entries as if
the latter were scalars. By a standard convention on linear operators, Θ ⊗ IH in (5) will be identified with the matrix
Θ.
For example, in the case when the system variables consist of the pairs of conjugate quantum mechanical position
qk and momentum pk := −i∂qk operators (with the Planck constant being appropriately normalised) [42, 55], with
1 6 k 6 n2 , they are defined on the Schwartz space [58] of rapidly decaying complex-valued functions on Rn/2 which
is dense in the Hilbert space L2 (Rn/2 ) of square integrable functions (and is invariant under those operators). In this
case, [q j , pk ] = iδ jk for all j, k = 1, . . . , 2n , where δ jk is the Kronecker delta, and the CCR matrix Θ of the system
variables q1 , . . . , qn/2 , p1 , . . . , pn/2 takes the form


1
1
0 In/2
Θ=
= J ⊗ In/2 ,
(6)
2 −In/2 0
2
where ⊗ is the Kronecker product of matrices, Ir denotes the identity matrix of order r, and use is made of the matrix


0 1
J :=
(7)
−1 0
which spans the space of antisymmetric matrices of order 2. In fact, any CCR matrix can be brought to the form (6)
by an appropriate linear transformation of the system variables.
The system under consideration is organised as a linear quantum stochastic system which models an open quantum
harmonic oscillator (OQHO) [13, 15] with n2 modes which interacts with an m-channel external bosonic field, where m

4

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

is even. The energetics of the system itself and its interaction with the fields is specified by the system Hamiltonian H
and the system-field coupling operators L1 , . . . , Lm which are, respectively, quadratic and linear functions of the system
variables:
 
L1
n
1
1 T
 .. 
L :=  .  = MX.
(8)
H := ∑ r jk X j Xk = X RX,
2 j,k=1
2
Lm
Here, R := (r jk )16 j,k6n ∈ Rn×n and M ∈ Rm×n are given matrices, with R being symmetric, which will be referred
to as the energy and coupling matrices. The evolution of the system variables in time t > 0 is governed by a linear
Hudson-Parthasarathy QSDE [23, 47]
dX = G (X)dt − i[X, LT ]dW = AXdt + BdW

(9)

whose structure (including its linearity in the case being considered) are clarified below. This QSDE is driven by the
vector
 
W1
 .. 
W :=  . 
Wm
of quantum Wiener processes W1 , . . . ,Wm which are time-varying self-adjoint operators on a symmetric Fock space F
[47, 49]. These operators represent the external fields and have a complex positive semi-definite Hermitian Ito matrix
Ω := (ω jk )16 j,k6m :
dW dW T = Ω dt,
Ω := Im + iJ.
(10)
Here,

J :=


0 Im/2
= J ⊗ Im/2
−Im/2 0

(11)

is an orthogonal real antisymmetric matrix of order m (so that J 2 = −Im ), which specifies CCRs for the quantum
Wiener processes as [dW, dW T ] = 2iJdt, and use is made of the matrix J from (7).
The Hilbert space H, which provides a common domain for the action of the system and external field variables,
has the tensor-product structure
H := H0 ⊗ F,
(12)
where H0 is the initial system space on which the system variables X1 (0), . . . , Xn (0) are defined. The matrices A ∈ Rn×n
and B ∈ Rn×m in (9) are expressed as
1
A := 2Θ (R + M T JM) = 2Θ R − BJBTΘ −1 ,
2

B := 2Θ M T

(13)

in terms of the energy and coupling matrices R and M from (8) and satisfy the physical realizability (PR) condition
[30, 56]:
AΘ +Θ AT + BJBT = 0.

(14)

Also, G in (9) denotes the Gorini-Kossakowski-Sudarshan-Lindblad generator [17, 36] (see also [1]), whose role is
similar to that of the infinitesimal generators of classical Markov diffusion processes [32]. More precisely, G is a
linear superoperator which specifies the drift term in the evolution of a system operator ξ (a function of the system
variables X1 , . . . , Xn ) as
dξ = G (ξ )dt − [ξ , hT ]dW.
(15)
Furthermore, G takes into account both the internal dynamics of the system and its interaction with the external fields
by acting on ξ as
G (ξ ) := i[H, ξ ] + D(ξ ),
(16)
where H is the system Hamiltonian in (8). The last term in (16) is associated with the system-field interaction and
involves the decoherence superoperator D which acts as
D(ξ ) :=

1
2

m

∑
j,k=1

ω jk (L j [ξ , Lk ] + [L j , ξ ]Lk ) .

(17)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

5

Here, ω jk are the entries of the quantum Ito matrix Ω in (10), and L1 , . . . , Lm are the coupling operators from (8). In
particular, by substituting L from (8) into (17) and using the CCRs (5), the entrywise application of D to the vector X
of the system variables in (3) leads to
D(X) :=

1
2

m

ω jk (L j [X, Lk ] + [L j , X]Lk )

∑
j,k=1

1
= 2Θ M T JMX = − BJBTΘ −1 X,
2

(18)

which gives rise to the field-related term 2Θ M T JM = − 12 BJBTΘ −1 in the matrix A in (13). The other part i[H, X] =
2Θ RX of the drift vector AX in the QSDE (9) comes from the quadratic nature of the Hamiltonian H in (8) in combination with the CCRs (5) and describes the internal dynamics of the system variables which they would have if the
system were isolated from the environment. Also, the representation of the dispersion matrix B = −i[X, LT ] of the
QSDE (9) in (13) follows from the linear dependence of the coupling operators in (8) on the system variables and the
CCRs (5).
In contrast to classical SDEs, the specific structure of the drift and diffusion terms of the QSDE (15) (and its
particular case (9)), described by (16)–(18), comes from the evolution
ξ (t) = U(t)† (ξ (0) ⊗ IF )U(t)

(19)

of the system operator ξ , with ξ (0) acting on the initial system space H0 . Here, U(t) is a time-varying unitary operator
which acts on the system-field space H in (12) and satisfies another QSDE


1
(20)
dU(t) = −U(t) i(H(t)dt + L(t)T dW (t)) + L(t)T Ω L(t)dt ,
2
with U(0) = IH . In accordance with the stochastic flow (19), the current Hamiltonian and coupling operators in (20)
are given by H(t) = U(t)† (H(0) ⊗ IF )U(t) and L(t) = U(t)† (L(0) ⊗ IF )U(t) and retain their quadratic-linear dependence (8) on the system variables X1 (t), . . . , Xn (t). This property follows from the fact that the map ζ 7→ U(t)† ζU(t),
which acts on operators ζ on H (and applies to vectors of such operators entrywise), is a unitary similarity transformation.
In view of (20), the unitary operator U(t) reflects an accumulated effect from internal driving forces of the system
and its interaction with the external fields over the time interval [0,t] and is adapted in the sense that it acts effectively
on the subspace H0 ⊗ Ft , where {Ft : t > 0} is the Fock space filtration. The QSDE (15) follows from (19) and (20)
due to the quantum Ito formula [23, 47] combined with (10), unitarity of U and commutativity between the forward Ito
increments dW (t) and adapted processes (including U) considered at time s 6 t. This commutativity is a consequence
of the continuous tensor-product structure [49] of the Fock space F.
The QSDEs (9), (15) and (20) correspond to a particular yet important scenario of quantum stochastic dynamics,
where there is no photon exchange between the fields, and the scattering matrix is an identity matrix, which effectively
eliminates the gauge processes from consideration. The presence of these processes and more general scattering matrices [23, 47] affects the dynamics of the unitary operator U and is taken into account in QSDEs which are used in the
quantum feedback network theory [18, 29].
Due to the linearity of the QSDE (9), the OQHO is employed extensively as a basic model in linear quantum
filtering and control [30, 46, 53] with quadratic performance criteria. If the initial system variables have finite second
moments, that is,
n


E X(0)T X(0) = ∑ E Xk (0)2 < +∞,
(21)
k=1

then, similarly to the classical case, this property is preserved in time by the linear dynamics (9) and leads to the finite
limit values of the first and second moments:

(22)
lim EX(t) = 0,
lim E X(t)X(t)T = P + iΘ ,
t→+∞

t→+∞

provided the matrix A in (13) is Hurwitz. Here, Eξ := Tr(ρξ ) is the expectation of a quantum variable ξ over the
density operator
ρ := ϖ ⊗ υ,
(23)
where ϖ is the initial system state on H0 , and υ is the vacuum state of the input bosonic fields on the Fock space F.
The matrix P in (22) is the infinite-horizon controllability Gramian
Z +∞

P :=
0

T

etA BBT etA dt

(24)

6

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

of the matrix pair (A, B), which is a unique solution of the following algebraic Lyapunov equation (ALE) due to the
matrix A being Hurwitz:
AP + PAT + BBT = 0.
(25)
The first two moments of the system variables, mentioned above, are part of the averaged behaviour of the system at a
particular instant and do not provide information on multi-point quantum correlations at different moments of time.
3 Multi-point Gaussian quantum states
We will first revisit the statistical properties of the system variables X1 (t), . . . , Xn (t) at one point in time t > 0. These
properties are encoded by the one-point quasi-characteristic function (QCF) Φ : R+ × Rn → C defined by (see, for
example, [10])
Φ(t, u) := EWu (t)
(26)
in terms of the Weyl operator (2) associated with X(t), with the averaging over the system-field state (23). Also, for
any t > 0, we denote by
Z
Z
t

Σ (t) :=

t

T

e(t−s)A BBT e(t−s)A ds =

0

T

esA BBT esA ds

(27)

0

the finite-horizon (over the time interval [0,t]) controllability Gramian of the pair (A, B), satisfying the Lyapunov ODE
Σ̇ = AΣ + Σ AT + BBT ,
with Σ (0) = 0. Since the matrix A is assumed to be Hurwitz, the function Σ is related to the matrix P in (24) by
T

Σ (t) = P − etA PetA → P,

as t → +∞.

(28)

The following lemma, which is given here for completeness, shows that the evolution of the QCF under the linear
QSDE is identical to that for classical Gaussian Markov diffusion processes (such as the Ornstein-Uhlenbeck process
[32]) produced by linear SDEs.
Lemma 1 Suppose the OQHO, governed by the linear QSDE (9), is driven by the input fields in the vacuum state.
Then the QCF Φ for the system variables in (26) satisfies a linear functional equation
T  − 1 kuk2
Φ(t, u) = Φ s, e(t−s)A u e 2 Σ (t−s)

(29)
√
n
T
for
√ any two moments of time t > s > 0 and all u ∈ R . Here, the function Σ is given by (27), and kvkK := v Kv =
| Kv| denotes a weighted Euclidean semi-norm of a vector v associated with a real positive semi-definite symmetric
matrix K.

Proof The solution of the linear QSDE (9) admits the decomposition
X(t) = e(t−s)A X(s) +Y (s,t)

(30)

for any t > s > 0, where Y (s,t) is an auxiliary vector of self-adjoint operators given by
Z t

Y (s,t) :=

e(t−τ)A BdW (τ).

(31)

s

Since the increments of the quantum Wiener process W commute with the current values of adapted processes, then
[X(s),Y (s,t)T ] =

Z t

T

[X(s), dW (τ)T ]BT e(t−τ)A = 0

(32)

s

for all t > s > 0. In fact, the entries of X(s) and Y (s,t) act on orthogonal subspaces of the system-field Hilbert space
H in (12). Furthermore, since the system and fields are in the product state (23), with the fields being in the vacuum
state, then X(s) and Y (s,t) are statistically independent. In combination with (30) and (32), this independence implies
that
T (e(t−s)A X(s)+Y (s,t))

Φ(t, u) = Eeiu

T e(t−s)A X(s)

= Eeiu

(t−s)AT

= Φ s, e

Ψ (s,t, u)

u Ψ (s,t, u),

(33)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

where
Ψ (s,t, u) := Eeiu

7

TY (s,t)

(34)

is the QCF for the vector Y (s,t). Now, in the vacuum state, the input fields have the quasi-characteristic functional
Eei

Rt

0

f (s)T dW (s)

1 R t | f (s)|2 ds
0

= e− 2

(35)

for any t > 0 and any locally square integrable function f : R+ → Rm . By substituting Y from (30) into (34) and using
(35), it follows that
Ψ (s,t, u) := Eeiu

T R t e(t−τ)A BdW (τ)
s

1 R t |BT e(t−τ)AT u|2 dτ
s

= e− 2

1 T R t−s eτA BBT eτAT dτu
0

= e− 2 u

− 12 kuk2Σ (t−s)

=e

,

(36)

where use is made of (27). Substitution of (36) into (33) leads to (29).

From (29), it follows that if the initial state ϖ of the system is Gaussian [48, 50] (that is, ln Φ(0, u) is a quadratic
function of u ∈ Rn ), then so is its reduced quantum state at subsequent moments of time t > 0. Furthermore, since
the matrix A is Hurwitz, the relations (28), (29) and the continuity limu→0 Φ(s, u) = Φ(s, 0) = 1 imply the pointwise
convergence of the QCF:


1
2
T  − 1 kuk2
u ∈ Rn ,
(37)
lim Φ(t, u) = lim Φ 0, etA u e 2 Σ (t) = e− 2 kukP ,
t→+∞

t→+∞

which holds regardless of whether the initial state is Gaussian (or whether (21) is satisfied). The relation (37) is
equivalent to the weak convergence [6, 10] of the reduced system state to a unique invariant Gaussian quantum state
with zero mean and the quantum covariance matrix P + iΘ specified by the matrix P in (24).
An infinite-dimensional extension of (29) is the quasi-characteristic functional of the quantum process X over the
time interval [0,t]:
e u) := Eei
Φ(t,

Rt

T
0 X(s) du(s)

Rt

T

sA

= Eei 0 (e X(0)+Y (0,s)) du(s)
 Zt T
 Rt Rs
T
= Φ 0, esA du(s) e− 0 ( 0 C(s,τ)du(τ)) du(s) ,

(38)

0

which is computed (similarly to that of a classical Gaussian random process) for any t > 0 and any function u : R+ →
Rn of locally bounded variation. The integrals in (38) are understood as appropriate (operator and vector-matrix)
versions of the Riemann-Stieltjes integral [68]. Also, use is made of the two-point covariance matrix
C(s, τ) := ReE(Y (0, s)Y (0, τ)T )
= e(s−τ)A Σ (τ) = C(τ, s)T ,

s > τ > 0,

(39)

for the quantum process Y (0, s) = 0s e(s−τ)A BdW (τ) from (30) and (31), where the function Σ is given by (27). The
covariance function C in (39) is the kernel function of the quadratic form
R

 Z t
2  Z
E
Y (0, s)T du(s)
=
0

du(s)TC(s, τ)du(τ)

[0,t]2
Z t Z s

=2

0

T
C(s, τ)du(τ) du(s).

0

In application to a piece-wise constant function u(t) = ∑16k6N: tk 6t vk with increments v1 , . . . , vN ∈ Rn at an increasing
e N , u) in (38) becomes
sequence of moments of time 0 6 t1 6 . . . 6 tN (with an arbitrary N = 1, 2, 3, . . .), the quantity Φ(t
the N-point QCF of the system variables:
N

T

ΦN (t1 , . . . ,tN ; v1 , . . . , vN ) := Eei ∑k=1 vk X(tk )
 N
 1 N
T
vTC(t ,t )v
−
= Φ 0, ∑ etk A vk e 2 ∑ j,k=1 j j k k .
k=1

(40)

8

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

Although (40) is a particular case of the quasi-characteristic functional (38), it is also possible to obtain (38) from (40)
by approximating the integrals with appropriate Riemann-Stieltjes sums. The multi-point QCFs ΦN with different
values of N are related to each other by
ΦN (t1 , . . . ,tN−1 ,tN−1 ; v1 , . . . , vN ) = ΦN−1 (t1 , . . . ,tN−1 ; v1 , . . . , vN−1 , vN−1 + vN ).

(41)

Furthermore, a multi-point extension of (29) is the recurrence relation
(tN −tN−1 )A

N−1 T

T

N−1 T

T (tN −tN−1 )A X(t
N−1 )

ΦN (t1 , . . . ,tN ; v1 , . . . , vN ) = Eei ∑k=1 vk X(tk )+ivN (e

X(tN−1 )+Y (tN−1 ,tN ))

= Eei ∑k=1 vk X(tk )+ivN e

T

EeivN Y (tN−1 ,tN )

 − 1 kvN k2Σ (t −t )
T
N N−1
= ΦN−1 t1 , . . . ,tN−1 ; v1 , . . . , vN−1 + e(tN −tN−1 )A vN e 2

(42)

which holds for all N = 2, 3, 4, . . ., where use is made of (27), (30) and (41). It follows from (42) by induction that
N ;v1 ,...,vN
) is a quadratic function of the vectors v1 , . . . , vN ∈ Rn whose coefficients depend on the time arguln ΦN (t1 ,...,t
t AT
N
Φ 0,∑k=1 e k

vk

ments t1 , . . . ,tN only through their differences tk − tk−1 . Therefore, since the matrix A is assumed to be Hurwitz, the
limit
−1 N
vTV (t −t )v
(43)
lim ΦN (t1 + t, . . . ,tN + t; v1 , . . . , vN ) = e 2 ∑ j,k=1 j j k k
t→+∞

is the QCF of a multi-point Gaussian quantum state with zero mean and quantum covariance matrix whose real part
is identical to the covariance matrix of a homogeneous Gaussian Markov diffusion process in Rn considered at the
moments of time 0,t2 −t1 , . . . ,tN −t1 . The covariance function V of this auxiliary classical process in (43) is computed
as
V (τ) = V (−τ)T = eτA P,
τ > 0,
(44)
where P is the matrix from (24). Due to the quantum nature of the setting under consideration, the Cn×n -valued
function
S(τ) := V (τ) + iΛ (τ) = S(−τ)∗
(45)
is a positive semi-definite Hermitian kernel (in the sense that so is the matrix (S(τ j − τk ))16 j,k6N for all τ1 , . . . , τN > 0
and N = 1, 2, 3, . . .). Here, (·)∗ := ((·))T denotes the complex conjugate transpose, and
Λ (τ) = −Λ (−τ)T = eτAΘ ,

τ > 0,

(46)

s,t > 0,

(47)

describes the two-point commutator matrix for the system variables:
[X(s), X(t)T ] = 2iΛ (s − t),

where use is made of (30)–(32). At the same time, S in (45) describes the covariance function of a classical stationary
Gaussian Markov diffusion process in the complex Euclidean space Cn , which will be discussed in Section 8.
As in the one-point case (37), the convergence of the multi-point QCFs in (43) holds regardless of whether the
initial system state is Gaussian (or whether (21) is satisfied). However, a Gaussian (but not necessarily invariant)
initial state gives rise to a multi-point Gaussian state associated with X(t1 ), . . . , X(tN ) at any future moments of time
tN > . . . > t1 > 0 for any N = 1, 2, 3, . . .. The quantum covariance matrices for such states are specified by the two-point
quantum covariance function of the system variables given by
cov(X(t), X(s)) := E(X(t)X(s)T ) − EX(t)EX(s)T
= e(t−s)A cov(X(s)),

(48)

cov(X(s)) := cov(X(s), X(s))
T

= esA cov(X(0))esA + Σ (s) + iΘ .

(49)

The relations (48) and (49) follow from the decomposition (30) and the fact that the entries of X(s) commute with and
are statistically independent of (and hence, uncorrelated with) those of Y (s,t) in (31).
While the mean square performance criteria in the quantum LQG control and filtering problems [13, 46, 61, 62]
are based on one-point second-order moments of the system variables, the multi-point Gaussian states involve the
covariance functions (48). These two-point moments are also involved in the risk-sensitive [11, 27, 28] cost functionals
which penalize higher-order moments of integral quantities.

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

9

4 Quadratic-exponential functional
Consider an adapted quantum process ϕ defined for any time t > 0 as the integral of a quadratic function ψ of the
system variables over the time interval [0,t]:
Z t

ϕ(t) :=

ψ(s)ds,

ψ(s) := X(s)T Π X(s).

(50)

0

Here, Π is a given real symmetric matrix of order n, so that ϕ(t) and ψ(t) are self-adjoint operators on the system-field
space H in (12). Furthermore, if Π < 0, then
n

ψ = ZTZ =

∑ Zk2 ,
k=1

 
Z1
 ..  √
Z :=  .  := Π X,

(51)

Zn

and both ϕ(t) and ψ(t) are positive semi-definite. Now, consider the exponential moment of the operator ϕ(t) from
(50) given by
Rt

Ξθ (t) := Eeθ ϕ(t) = Ee

T
0 X(s) Π X(s)ds

,

(52)

where θ is a nonnegative real-valued parameter (which is assumed to be sufficiently small in order for Ξθ (t) to be
finite). The dependence of Ξθ (t) on the matrix Π through (50) is omitted for brevity. The quadratic-exponential
functional (QEF) in (52) can be regarded as an alternative to the quantum risk-sensitive performance criterion in [11,
27, 28]. The latter was defined as a weighted mean square value of a time-ordered exponential satisfying an operator
differential equation (see [28, Eqs. (19)–(21)]). In fact, the QEF is a straightforward quantum version of its classical
predecessors [5, 26, 64]. For any given time t > 0, the quantity Ξθ (t) is the moment-generating function for the quantum
variable ϕ(t) in the sense that
E(ϕ(t)k ) = ∂θk Ξθ (t)

(53)

θ =0

for any positive integer k for which the moment exists. In particular, if the system variables have finite second moments
(21), then the asymptotic behaviour of the QEF for small values of θ is described by
Z t

Ξθ (t) = 1 + θ

hΠ , P(s)ids + o(θ ),

θ → 0,

(54)

0

where hK, Ni := Tr(K ∗ N) denotes the Frobenius inner product [22] of real or complex matrices K and N. Here,
P(t) := ReE X(t)X(t)T



(55)

is the real part of the matrix E(XX T ) = P + iΘ of second moments of the system variables governed by the Lyapunov
ODE
Ṗ = AP + PAT + BBT

(56)

whose initial condition satisfies the generalized Heisenberg uncertainty principle P(0) + iΘ < 0 (see, for example,
[21]). Moreover, the first two terms on the right-hand side of (54) provide a lower bound for the QEF. Indeed, for
any fixed but otherwise arbitrary time t > 0, the self-adjoint operator ϕ(t) can be regarded as a classical real-valued
random variable with a probability distribution Et .1 Hence,
Z +∞

Ξθ (t) =
>

eθ x Et (dx)

−∞
Z +∞
−∞

(1 + θ x)Et (dx) = 1 + θ Eϕ(t)
Z t

= 1+θ

Z t

Eψ(s)ds = 1 + θ
0

hΠ , P(s)ids,

(57)

0

1 This distribution is related by E (A) := EP (A) to the spectral measure P of ϕ(t), which is a projection-valued measure on the σ -algebra B
t
t
t
T
of Borel subsets of the real line satisfying Pt (A)Pt (B) = Pt (A B) for all A, B ∈ B and the resolution of the identity property Pt (R) = I ; see, for
example, [21].

10

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

where use is also made of the inequality ev > 1 + v on the real line, and P is the function given by (55) and (56).
Furthermore, in combination with the first equality in (57) and convexity of the exponential function, the Jensen
inequality leads to
Z +∞

eθ x Et (dx)

ln Ξθ (t) = ln

−∞
Z +∞

>θ

−∞
Z t

=θ

xEt (dx) = θ Eϕ(t)
Z t

Eψ(s)ds = θ
0

hΠ , P(s)ids.

(58)

0

In view of (57) and (58), the QEF Ξθ (t) in (52) and its logarithm provide upper bounds for the mean square cost
functional which is used in the quantum LQG control and filtering problems [13, 30, 46]. An appropriate scaling
represents the QEF as a convex combination of the moments from (53):
e−θ Ξθ (t) = e−θ

+∞

θk
E(ϕ(t)k ).
k=0 k!

(59)

∑

k

Its coefficients θk! e−θ constitute the Poisson probability mass function [57] with intensity parameter θ > 0 and achieve
their maximum at k = bθ c, where b·c is the floor function. Therefore, θ determines the order of the moment of ϕ(t),
which is endowed with the largest weight in (59) and hence, is penalized most by the QEF as a cost functional. In
comparison with the quadratic cost functionals in LQG control and filtering problems, the QEF leads to minimizing
not only the one-point second-order moments of the system variables but also their multi-point higher-order moments,
with the latter being subject to a stronger penalty for larger values of θ . A similar consideration applies to the cumulantgenerating function
+∞

θr
Kr (ϕ(t))
r=1 r!


θ
= θ Eϕ(t) + var(ϕ(t)) + O(θ 3 ),
2

ln Ξθ (t) =

∑

as θ → 0,

(60)

for the self-adjoint quantum variable ϕ(t) as a classical random variable (with the distribution Et at a given moment of
time t). Here, the variance var(φ ) := E(φ 2 ) − (Eφ )2 = K2 (φ ) of a random variable φ is the second of the cumulants
Kr (φ ) := ∂vr ln Eevφ

v=0

= Pr (Eφ , . . . , E(φ r )).

(61)

For any r > 1, the rth cumulant Kr (φ ) is expressed in terms of the first r moments Eφ , . . . , E(φ r ) through the polynomial
r

(−1)k
k
k=1

Pr (µ1 , . . . , µr ) = −r! ∑

k

∑

j1 ,..., jk >1: j1 +...+ jk =r

µ js
s=1 js !

∏

r

(−1)k
k
k=2

= µr − r! ∑

k

∑

j1 ,..., jk >1: j1 +...+ jk =r

µ js
,
s=1 js !

∏

(62)

whose coefficients do not depend on the probability distribution of φ . In particular, the first three polynomials are
given by
P1 (µ1 ) = µ1 ,
P2 (µ1 , µ2 ) = µ2 − µ12 ,
P3 (µ1 , µ2 , µ3 ) = µ3 − 3µ1 µ2 + 2µ13 .
In view of (50), the first two terms on the right-hand side of (60) involve the second and fourth moments of the system
variables over the time interval [0,t] and will be employed in Section 5 for approximating the QEF.
As mentioned above, the QEF (52) differs from its predecessor in [27] and [28, Eqs. (19)–(21)] which was used for
measurement-based risk-sensitive quantum control problems and defined in terms of an operator differential equation.
Nevertheless, the QEF satisfies an integro-differential equation with somewhat similar features, which will be provided

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

11

by Theorem 1 below. For its formulation, we will need an auxiliary integral operator which is introduced by the
following lemma. To this end, for any time t > 0, let ft denote a linear map which acts as
Z t


ft (α, β ) :=Re X(t)T α(σ )X(σ )dσ
0

Z

+

[0,t]2

X(σ )T β (σ , τ)X(τ)dσ dτ

(63)

on a pair of locally integrable matrix-valued functions α : R+ → Rn×n and β : R2+ → Rn×n , the second of which
satisfies
β (σ , τ)T = β (τ, σ )

(64)

for all σ , τ > 0. Here, the real part Re(·) is extended to operators as Reξ := 12 (ξ + ξ † ). In view of (64), the image
ft (α, β ) in (63) is a self-adjoint quantum variable on H which depends in a quadratic fashion on the past history of
the system variables of the OQHO over the time interval [0,t].
Lemma 2 For any given t > 0, the composition of the map ft in (63) with the commutator adϕ(t) (·) := [ϕ(t), ·],
associated with the process ϕ in (50), can be represented as
adϕ(t) ft = ift Γt .

(65)

e , βe) = Γt (α, β )
Here, Γt is a linear integral operator which maps the function pair (α, β ) in (63) and (64) to the pair (α
of such functions as
e (σ ) := − 4
α

Z t

α(s)Λ (s − σ )dsΠ ,

(66)

0

βe(σ , τ) :=2(ΠΛ (σ − t)α(τ) − α(σ )TΛ (t − τ)Π )
Z t


ΠΛ (σ − s)β (s, τ) − β (σ , s)Λ (s − τ)Π ds,

(67)

where Λ is the two-point commutator function of the system variables from (46) and (47).



+4
0

Proof Note that
1
1
[ξ , Reη] = ([ξ , η] + [ξ , η † ]) = ([ξ , η] + [η, ξ † ]† )
2
2
1
†
= ([ξ , η] − [ξ , η] ) = iIm[ξ , η]
2
for any self-adjoint operator ξ and an arbitrary operator η, where Imζ :=
imaginary part Im(·). A combination of (50), (63) and (68) implies that

1
2i (ζ

(68)

− ζ † ) is the operator extension of the

Z t

[ϕ(t), ft (α, β )] =

[ψ(s), ft (α, β )]ds
Z t
i

=
ψ(s), Re X(t)T α(σ )X(σ )dσ ds
0
0
Z
Z th
i
+
ψ(s),
X(σ )T β (σ , τ)X(τ)dσ dτ ds
0

Z th

[0,t]2

0

Z

=iIm

[0,t]2

Z

+

[0,t]3

[ψ(s), X(t)T α(σ )X(σ )]dsdσ

[ψ(s), X(σ )T β (σ , τ)X(τ)]dsdσ dτ.

(69)

Now, the quadratic polynomials ∑rj=1 ∑sk=1 γ jk X j (σ j )Xk (τk ) of the system variables of the OQHO, considered at arbitrary moments of time σ1 , . . . , σr and τ1 , . . . , τs , with complex coefficients γ jk (or integrals of such polynomials), form
a Lie algebra with respect to the commutator; see Appendix A and references therein. Furthermore, the two-point

12

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

CCRs (47) allow (A5) of Lemma 4 to be applied to the quadratic functions of the system variables on the right-hand
side of (69) as
[ψ(s), X(t)T α(σ )X(σ )] = [X(s)T Π X(s), X(t)T α(σ )X(σ )]

T 


0
ΠΛ (s − t)α(σ ) X(s)
X(s)
= 4i
X(t)
X(σ )
α(σ )Λ (s − σ )T Π
0

= 4i X(s)T ΠΛ (s − t)α(σ )X(σ ) − X(t)T α(σ )Λ (σ − s)Π X(s) ,
T

T

(70)

T

[ψ(s), X(σ ) β (σ , τ)X(τ)] = [X(s) Π X(s), X(σ ) β (σ , τ)X(τ)]

T 


0
ΠΛ (s − σ )β (σ , τ) X(s)
X(s)
= 4i
X(σ )
X(τ)
β (σ , τ)Λ (s − τ)T Π
0

= 4i X(s)T ΠΛ (s − σ )β (σ , τ)X(τ) − X(σ )T β (σ , τ)Λ (τ − s)Π X(s) ,

(71)

where the first equality from (46) has also been taken into account. Substitution of (70) and (71) into (69) leads to
Z

[ϕ(t), ft (α, β )] =4iRe

[0,t]2

Z

+ 4i

[0,t]3


X(s)T ΠΛ (s − t)α(σ )X(σ ) − X(t)T α(σ )Λ (σ − s)Π X(s) dsdσ

X(s)T ΠΛ (s − σ )β (σ , τ)X(τ) − X(σ )T β (σ , τ)Λ (τ − s)Π X(s) dsdσ dτ

Z

= − 4iRe X(t)T

[0,t]2

Z

+ 4iRe

[0,t]2

Z

+ 4i

[0,t]3

α(s)Λ (s − σ )Π X(σ )dsdσ



X(σ )T ΠΛ (σ − t)α(τ)X(τ)dσ dτ


X(σ )T ΠΛ (σ − s)β (s, τ) − β (σ , s)Λ (s − τ)Π X(τ)dsdσ dτ,

(72)

where use is also made of appropriate permutations of the integration variables. By grouping the terms on the righthand side of (72) and comparing the result with (63), it follows that
e , βe) = ift (Γt (α, β )),
[ϕ(t), ft (α, β )] = ift (α

(73)

e and βe are related to α and β through the linear operator Γt specified by (66) and (67). The
where the functions α
e
function β in (67) inherits the property (64) due to the symmetry of the matrix Π and the first equality in (46). The
fulfillment of (73) for arbitrary functions α and β , described above, establishes (65).

Lemma 2 shows that the class of quantum variables ft (α, β ), described by (63) for admissible functions α and β
is invariant under the action of iadϕ(t) . Moreover, (65) extends as
f (iadϕ(t) )ft = ft f (−Γt )

(74)

to entire functions f of a complex variable evaluated at operators. The nature of the operator f (−Γt ) in (74) is simpler
than that of f (iadϕ(t) ). Indeed, Γt (for any fixed but otherwise arbitrary t > 0) in (66), (67) is a bounded operator on
the Banach space of pairs of integrable Rn×n -valued functions on the time interval [0,t], described above, and hence,
so also is the operator f (−Γt ). This boundedness follows from the uniform boundedness of the function Λ in (46) due
to the matrix A being Hurwitz.
Also, for the purposes of the following theorem, we will need a time-varying change of the density operator in
(23):
θ
θ
1
e 2 ϕ(t) ρe 2 ϕ(t) .
(75)
ρθ ,t :=
Ξθ (t)
The property that ρθ ,t is a density operator is ensured by its self-adjointness, positive semi-definiteness and the unit
trace

 θ
θ
1
Trρθ ,t =
Tr e 2 ϕ(t) ρe 2 ϕ(t)
Ξθ (t)

1
1
=
Tr ρeθ ϕ(t) =
Eeθ ϕ(t) = 1
Ξθ (t)
Ξθ (t)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

13

in view of (52). The expectation over the modified density operator ρθ ,t in (75) is denoted by Eθ ,t and computed as
Eθ ,t ξ := Tr(ρθ ,t ξ )

θ
θ
1
=
Tr e 2 ϕ(t) ρe 2 ϕ(t) ξ
Ξθ (t)

θ
θ
1
=
E e 2 ϕ(t) ξ e 2 ϕ(t) .
Ξθ (t)

(76)

In the case θ = 0, the definitions (75) and (76) reproduce the original density operator ρ0,t = ρ and the original
expectation E0,t = E.
Theorem 1 As a function of time, the QEF Ξθ (t) in (52), associated with the OQHO (9), satisfies the integrodifferential equation

∂t ln Ξθ (t) = θ Eθ ,t ψ(t) +ϒθ (t) ,
(77)
where Eθ ,t is the modified expectation (76), and the process ψ is given by (50). Here, ϒθ (t) is a time-varying selfadjoint operator which depends in a quadratic fashion on the past history of the system variables according to (63)
as
θ
ϒθ (t) := ft (αθ ,t , βθ ,t )
2
Z t


θ
=
Re X(t)T αθ ,t (σ )X(σ )dσ
2
0
Z

+



T

[0,t]2

X(σ ) βθ ,t (σ , τ)X(τ)dσ dτ .

(78)

Also, the pair of Rn×n -valued functions αθ ,t and βθ ,t (the second of which satisfies the symmetry property (64)) is
computed as the image
θ 
bt , 0)
Γt (α
(79)
(αθ ,t , βθ ,t ) = f
2

bt , 0) under the linear operator f θ2 Γt , where the integral operator Γt is described in Lemma 2. Also, the
of the pair (α
bt is given by
function α
bt (σ ) := −8ΠΛ (t − σ )Π
α
(80)
in terms of the two-point commutator function Λ of the system variables from (46) and (47). Furthermore,

f (z) :=

0
sin z−z
z2

+∞
if z = 0
z2k−1
= ∑ (−1)k
otherwise k=1
(2k + 1)!

(81)

is an entire odd function of a complex variable which is evaluated in (79) at a bounded operator θ2 Γt .



Proof Application of the Magnus lemma on the differentiation of operator exponentials [20, 33, 37, 40, 65] leads to


(eθ ϕ(t) ) = θ

Z 1

eλ θ ϕ(t) ϕ̇(t)e(1−λ )θ ϕ(t) dλ

0
θ

= θ e 2 ϕ(t)
= θe

Z

1
2

− 12
Z 1
θ ϕ(t)
2
2
− 12

θ

eλ θ ϕ(t) ψ(t)e−λ θ ϕ(t) dλ e 2 ϕ(t)
θ

eλ θ adϕ(t) (ψ(t))dλ e 2 ϕ(t)

θ

θ
θ
= θ e 2 ϕ(t) sinhc adϕ(t) (ψ(t))e 2 ϕ(t) .
2
Here, ϕ̇(t) = ψ(t) is the time derivative of the process ϕ in (50), and
eµϕ ψe−µϕ = eµadϕ (ψ) =

+∞

µk k
adϕ (ψ),
k=0 k!

∑

(82)

14

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James
k times

z
}|
{
where
:= [ϕ, [ϕ, . . . [ϕ, · ] . . .] is the k-fold iterate of the commutator adϕ(t) (·) (with ad0ϕ being the identity map
by the standard convention). Also,

Z
1 1 λz
1 if z = 0
e dλ = sinh z
sinhc(z) :=
2 −1
z otherwise
adkϕ (·)

+∞

=

z2k

(83)

∑ (2k + 1)!

k=0

is a hyperbolic version of the sinc function (in the sense that sinhc(iz) = sinc(z)). This is an entire even function of a
complex variable z, which is evaluated in (82) at the superoperator θ2 adϕ . Since sinhc(0) = 1 in view of (83), then

θ
sinhc adϕ(t) (ψ(t)) = ψ(t) +ϒθ (t),
2

(84)

where
+∞

ϒθ (t) :=

(θ /2)2k

∑ (2k + 1)! ad2kϕ(t) (ψ(t))

k=1


θ θ
= g adϕ(t) ([ϕ(t), ψ(t)]).
2 2

(85)

Here,

g(z) :=

0
sinhc(z)−1
z

+∞
if z = 0
z2k−1
=∑
otherwise k=1 (2k + 1)!

(86)

is an entire odd function on the complex plane, which is applied to the superoperator θ2 adϕ(t) in (85), with the resulting
superoperator acting on the (skew-Hermitian) quantum variable [ϕ(t), ψ(t)]. The self-adjointness of the operator ϒθ (t)
follows directly from (85). Indeed, since adiξ (η) = i[ξ , η] is a self-adjoint operator for any self-adjoint operators ξ
and η, then so also is ad2ξ (η) = −ad2iξ (η), and, by induction, ad2k
ξ (η) is self-adjoint for all k = 1, 2, 3, . . .. Therefore,
the series (85) consists of self-adjoint operators due to the self-adjointness of ϕ(t) and ψ(t) in (50). Now, substitution
of (84) into (82) yields


θ

θ

(eθ ϕ(t) ) = θ e 2 ϕ(t) (ψ(t) +ϒθ (t))e 2 ϕ(t) .

(87)

By taking the expectation on both sides of (87) and using (52) together with the modified expectation from (76), it
follows that

∂t Ξθ (t) = E (eθ ϕ(t) )

 θ
θ
= θ E e 2 ϕ(t) (ψ(t) +ϒθ (t))e 2 ϕ(t)

= θ Ξθ (t)Eθ ,t ψ(t) +ϒθ (t) .
(88)
Division of both parts of (88) by the QEF Ξθ (t) leads to its logarithmic derivative in (77). We will now compute the
process ϒθ in (85). To this end, use will be made of the Lie-algebraic property of quadratic polynomials of system
variables of the OQHO, which was already employed in the proof of Lemma 2. More precisely, in view of the twopoint CCRs (47), application of (A5) from Lemma 4 to the quadratic functions of the system variables in (50) leads
to
[ψ(s), ψ(t)] = [X(s)T Π X(s), X(t)T Π X(t)]

T 


X(s)
0
γ(s − t) X(s)
= 4i
X(t)
−γ(t − s)
0
X(t)

= 4i X(s)T γ(s − t)X(t) − X(t)T γ(t − s)X(s)

= −8iRe X(t)T γ(t − s)X(s)

(89)

for any s,t > 0. Here,
γ(τ) := ΠΛ (τ)Π = −γ(−τ)T ,

(90)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

15

and use is also made of (46) in combination with the symmetry of the matrix Π and the identity (ξ T Nη)† = ηN T ξ
which holds for vectors ξ and η of self-adjoint operators and an appropriately dimensioned real matrix N. In view of
(50), it follows from (89) that
Z t

[ϕ(t), ψ(t)] =

[ψ(s), ψ(t)]ds
Z t


= −8iRe X(t)T γ(t − s)X(s)ds
0

0

bt , 0),
= ift (α

(91)

bt given by (80) in view of (90).
where the last equality is obtained by comparison with (63), with the function α
Substitution of (91) into (85) leads to

θ θ
bt , 0))
ϒθ (t) = g adϕ(t) (ift (α
2 2




θ
θ
bt , 0)
= i ft g i Γt (α
2
2

θ  θ 
bt , 0) .
= ft f
Γt (α
(92)
2
2
Here, use is made of Lemma 2 (and its corollary (74)), with the function f in (81) being related to g in (86) by
f (z) := ig(iz). In view of the notation (79), the representation (92) establishes (78), thus completing the proof.

In the classical case (obtained in the limit as Θ → 0 and J → 0), the operators ϕ(t) and ψ(t) in (50) are real-valued
(and hence, commuting) functions of time, and the process ϒθ vanishes. In this case, (77) reduces to


θ
∂t ln Ξθ (t) = θ Eθ ,t ψ(t) =
E eθ ϕ(t) ψ(t)
Ξθ (t)
and corresponds to the derivative of the usual exponential function. Returning to the quantum setting, we note that
(87) resembles the operator differential equation [28, Eq. (19)]. However, unlike the latter, (87) produces a self-adjoint
solution, and its right-hand side involves an additional process ϒθ which depends on the past history of the system
variables. More precisely, the operator ϒθ (t) stores part of the quantum information on the past history of the system

variables that is relevant to the time derivative (eθ ϕ(t) ) .
5 A quartic approximation of the quadratic-exponential functional
As mentioned in the previous section in regard to the Taylor series expansion (60), the QEF can be approximated by
the quantity


θ
Fθ (t) := θ Eϕ(t) + var(ϕ(t))
(93)
2
which takes into account only the first four moments of the system variables over the time interval [0,t]. In view of
(57) and (58), the term Eϕ(t) in (93) is organised as a mean square cost functional in the coherent quantum LQG
control problems [46]. This quadratic term is revisited for completeness in the following theorem whose main result
is the computation of the second term in (93).
Theorem 2 Suppose the matrix A in (13) is Hurwitz, and the OQHO (9), driven by vacuum fields, is initialised at the
invariant Gaussian state. Then the mean value and variance of the quantum process ϕ in (50) can be computed as
Eϕ(t) = hΠ , Pit,
Z t

var(ϕ(t)) = 4

(94)

D
E
T
(t − τ) Π , eτA (PΠ P +Θ ΠΘ )eτA dτ

(95)

0

for any t > 0, where the matrix P is given by (24). The infinite-horizon asymptotic behavior of this variance is described
by
1

lim
var(ϕ(t)) = 4hΠ , T i = 4hQ, PΠ P +Θ ΠΘ i,
(96)
t→+∞ t
where the matrices T and Q are the unique solutions of the ALEs
AT + TAT + PΠ P +Θ ΠΘ = 0,
T

A Q + QA + Π = 0.

(97)
(98)


16

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

Proof Since the invariant Gaussian quantum state has zero mean and the covariance matrix P + iΘ , the relation (94)
follows from (50) and the equality
Eψ(t) = Π , E(X(t)X(t)T ) = hΠ , Pi
(99)
for any t > 0. Also, (50) allows the variance of ϕ(t) to be expressed in terms of the quantum covariance function of ψ
as
Z
var(ϕ(t)) =
cov(ψ(σ ), ψ(τ))dσ dτ.
(100)
[0,t]2

Application of Lemma 5 of Appendix B to the two-point Gaussian quantum state for the system variables at the
moments of time σ and τ leads to
cov(ψ(σ ), ψ(τ)) = cov(X(σ )T Π X(σ ), X(τ)T Π X(τ))
= 2 Π , S(σ − τ)Π S(σ − τ)T ,

(101)

where use is also made of the symmetry of the matrix Π . Here, in accordance with (44)–(49), the corresponding
quantum covariance matrix of the system variables
(
e(σ −τ)A (P + iΘ ) if σ > τ > 0
T
E(X(σ )X(τ) ) =
T
(P + iΘ )e(τ−σ )A if τ > σ > 0
= V (σ − τ) + iΛ (σ − τ) = S(σ − τ)

(102)

depends on the difference σ − τ since the system is assumed to be initialised at the invariant state. A combination of
(101) with an appropriate transformation of the integration variables in (100) leads to
Z

var(ϕ(t)) = 2

[0,t]2
Z t

=2
−t
Z t

Π , S(σ − τ)Π S(σ − τ)T dσ dτ

(t − |τ|) Π , S(τ)Π S(τ)T dτ

(t − τ)hΠ , S(τ)Π S(τ)T + S(−τ)Π S(−τ)T idτ

=2
0

Z t

=4

(t − τ) Π , Re S(τ)Π S(τ)T



dτ.

(103)

0

Here, the symmetry of the matrix Π has been used together with the property S(−τ) = S(τ)∗ (see also (44)–(46)),
whereby
E D
E
D
Π , S(−τ)Π S(−τ)T = Π , S(τ)∗ Π S(τ) = S(τ)Π , Π S(τ)
D
E
= Π S(τ), S(τ)Π = hΠ , S(τ)Π S(τ)T i
for any τ ∈ R. By substituting (102) into (103) and using the relation (P + iΘ )T = P − iΘ in view of the symmetry of
P and antisymmetry of Θ , it follows that
Z t

 T
(t − τ)hΠ , eτA Re (P + iΘ )Π (P − iΘ ) eτA idτ
0
Z t
D
E
T
= 4 (t − τ) Π , eτA (PΠ P +Θ ΠΘ )eτA dτ,

var(ϕ(t)) = 4

0

which establishes (95). It now remains to prove (96). To this end, the representation (103) implies that, for any t > 0,
1
var(ϕ(t)) = 4
t
where

Z +∞
0

χt (τ) Π , Re S(τ)Π S(τ)T



dτ,

(104)


τ
χt (τ) := max 0, 1 −
(105)
t
is a finite support function which is bounded by and converges to 1 as t → +∞ for any given τ > 0. Since (due to A
being Hurwitz) the covariance function S(τ) in (102) is square integrable over τ, application of Lebesgue’s dominated
convergence theorem to (104) leads to
Z +∞
1


lim
var(ϕ(t)) = 4
Π , Re S(τ)Π S(τ)T dτ = 4 hΠ , T i ,
(106)
t→+∞ t
0

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

17

thus proving the first equality in (96), where the matrix
Z +∞

T :=

Z

Re S(τ)Π S(τ)T dτ =

0

+∞

T

eτA (PΠ P +Θ ΠΘ )eτA dτ

(107)

0

satisfies the ALE (97). On the other hand, application of a duality argument allows the integral in (106) to be evaluated
as
Z +∞ D
Z +∞
E

T
Π , eτA (PΠ P +Θ ΠΘ )eτA dτ
Π , Re S(τ)Π S(τ)T dτ =
0
0
Z +∞ D
E
T
=
eτA Π eτA , PΠ P +Θ ΠΘ dτ
0

= hQ, PΠ P +Θ ΠΘ i ,
T

which establishes the second of the equalities in (96), with the matrix Q := 0+∞ eτA Π eτA dτ satisfying the ALE (98).

In the limiting case Θ = 0, when the system variables commute with each other, the relation (101) reduces to the
representation
cov(ξ T Π ξ , η T Π η) = 2 hΠ , cov(ξ , η)Π cov(η, ξ )i
(108)
R

for the covariance of quadratic forms of jointly Gaussian classical random vectors ξ and η with zero mean values (see,
for example, [41, Lemmas 2.3 and 6.2] and [59, Lemma 6]). The noncommutative quantum nature of the setting under
consideration with Θ 6= 0 manifests itself in (95)–(97) through the term Θ ΠΘ .
Also, in application to the quartic approximation (93) of the QEF, Theorem 2 yields the following infinite-horizon
asymptotic behaviour:

1
1
lim
Fθ (t) = hΠ , P + 2θ T i.
(109)
θ t→+∞ t
Here, T is the matrix from (107) which takes into account the two-point quantum statistical correlations (102) between
the system variables. The relation (109) shows that, in the infinite-horizon limit, the quartic term θ2 var(ϕ(t)) in (93)
can be neglected in comparison with the quadratic term Eϕ(t) only if the parameter θ is small enough in the sense
that
1 hΠ , Pi
θ  θ0 :=
.
(110)
2 hΠ , T i
6 Asymptotic growth rates of the higher-order cumulants
Similarly to the quadratic and quartic terms of the Taylor series expansion of the QEF in (60), the higher-order cumulants Kk (ϕ(t)) also grow asymptotically linearly with time t. For the formulation of the following theorems, we will
need an auxiliary transformation

N(τ) if k = 0
[k]
N (τ) :=
(111)
N(−τ)T if k = 1
of a matrix-valued function N (of a real variable), which is the identity transformation in the case k = 0 (that is,
N [0] = N). Also, we denote by
r
−
→
(112)
∏ Nk := N1 × . . . × Nr
k=1

the rightward-ordered product of operators or appropriately dimensioned matrices (with the order being essential
because of noncommutativity). Furthermore, let

0 if a 6 b
(113)
I(a, b) :=
1 if a > b
indicate an inversion in the pair of integers (a, b). The inversion indicator I extends to consecutive pairs in an r-tuple
of integers as
I(c1 , c2 , . . . , cr ) := (I(c1 , c2 ), I(c2 , c3 ), . . . , I(cr−1 , cr )).
(114)
Also, for any binary (r − 2)-index γ ∈ {0, 1}r−2 , we denote by ∆r,γ the number of permutations c := (c1 , . . . , cr−1 ) of
the set {1, . . . , r − 1} with the inversion indicator γ:
∆r,γ := #{c is a permutation of {1, . . . , r − 1} such that I(c) = γ}.

(115)

18

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

Since such permutations, considered for all possible γ ∈ {0, 1}r−1 , form a partitioning of the set of (r − 1)! permutations of {1, . . . , r − 1}, then
(116)
∑ ∆r,γ = (r − 1)!.
γ∈{0,1}r−2

Theorem 3 Suppose the matrix A in (13) is Hurwitz, and the OQHO (9), driven by vacuum fields, is initialised at the
invariant Gaussian state. Then for any r > 2 and any time t > 0, the rth cumulant (61) of the process ϕ in (50), with
Π < 0, can be represented as
Kr (ϕ(t)) = 2r−1

Z

∑

∆r,γ

γ∈{0,1}r−2

[0,t]r

r−1


−
→

Tr Π S(t1 − t2 ) ∏ Π S[γ j ] (t j − t j+1 ) Π S(t1 − tr )T dt1 × . . . × dtr ,

(117)

j=2

where S is the steady-state quantum covariance function of the system variables in (102), and use is made of the
notation (111) and (112). The sum in (117) is over binary (r − 2)-indices γ := (γ2 , . . . , γr−1 ) ∈ {0, 1}r−2 , and the
corresponding coefficients ∆r,γ are described by (114)–(116).

Proof In what follows, for any given r = 1, 2, 3, . . ., a permutation
κ := (κ1 , . . . , κ2r ) = ( j1 , k1 , . . . , jr , kr )

(118)

of the integers 1, . . . , 2r is called regular if it satisfies
j1 < j2 < . . . < jr−1 < jr ,

j1 < k1 , j2 < k2 , . . . , jr < kr .

(119)

The set of regular permutations consists of (2r)!
r!2r = (2r − 1)!! elements and is denoted by Pr . There is a one-to-one
correspondence between Pr and the set of all partitions of {1, . . . , 2r} into two-element subsets. Note that 1 is a fixed
point for every regular permutation since j1 = 1. For example, the set of regular permutations of {1, 2, 3, 4} is given
by
P2 = {(1, 3, 2, 4), (1, 4, 2, 3), (1, 2, 3, 4)}.
The Wick-Isserlis theorem2 [31, Theorem 1.28 on pp. 11–12] (see also [24] and [39, p. 122]) allows the product
moment of self-adjoint quantum variables ξ1 , . . . , ξ2r , which satisfy CCRs and are in a zero-mean Gaussian state, to
be expressed in terms of their covariances as
2r

−
→

r

E ∏ ξk =

∑ ∏ E(ξ js ξks ).

(120)

κ∈Pr s=1

k=1

By using (50) and (51) and applying (120) to the multi-point zero-mean Gaussian state of the entries of the process Z
at different moments of time, it follows that
 Z
E(ϕ(t)r ) = E

t n

∑ Z` (τ)2 dτ

r 

0 `=1
n

=
=
=

∑

r

Z

r
`1 ,...,`r =1 [0,t]
Z
n

∑

−
→
E ∏ (Z`k (tk )2 )dt1 × . . . × dtr
k=1

r

∑ ∏ K`d( js ) `d(ks ) (td( js ) − td(ks ) )dt1 × . . . × dtr

r
`1 ,...,`r =1 [0,t] κ∈Pr s=1
Z
n
r

∑

∑ ∏ K`d( js ) `d(ks ) (td( js ) − td(ks ) )dt1 × . . . × dtr .

r
κ∈Pr [0,t] `1 ,...,`r =1 s=1

(121)

Here, K denotes the steady-state quantum covariance function of Z in (51) which is related to that of the process X in
(102) by
√
√
E(Z(σ )Z(τ)T ) = Π S(σ − τ) Π
=: K(σ − τ) = (Kab (σ − τ))16a,b6n .
2

it is also used in Appendix B

(122)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

19

Also, the function
d(k) := dk/2e

(123)

in (121) maps the set {1, . . . , 2r} onto {1, . . . , r}, where d·e denotes the ceiling function. For any j = 1, . . . , r, the
inverse image of the singleton { j} under the map d is the two-element set d −1 ( j) = {2 j − 1, 2 j}. For any κ ∈ Pr
from (118), the 2r-tuple
ν := (ν1 , . . . , ν2r ) := d(κ)
:= (d( j1 ), d(k1 ), . . . , d( jr ), d(kr ))

(124)

describes a permutation of the multi-set {1, 1, 2, 2, . . . , r − 1, r − 1, r, r}. Since the function d in (123) is nondecreasing,
then ν inherits from κ the property (119) in a nonstrict form:
ν1 6 ν3 6 . . . 6 ν2r−3 6 ν2r−1 ,

ν1 6 ν2 , ν3 6 ν4 , . . . , ν2r−1 6 ν2r .

(125)

Since the map d is not injective, some elements κ ∈ Pr are mapped to the same permutation ν in (124), so that
the inverse image d −1 (ν) := {κ ∈ Pr : d(κ) = ν} of ν ∈ d(Pr ) consists, in general, of more than one regular
permutation. In fact, only the identity permutation has a unique preimage: d −1 (1, 1, 2, 2, . . . , r, r) = (1, 2, . . . , 2r −1, 2r).
For example, the case r = 3 is illustrated by Tabs. 1 and 2 which also provide the numbers #d −1 (ν) of elements in the
inverse images. For any r > 1, the relation (121) can be represented in terms of the permutations ν in (124) (with their
Table 1 The set P3 of regular permutations κ of {1, . . . , 6} and their images ν := d(κ) under the map d in (123)
No.

κ1

κ2

κ3

κ4

κ5

κ6

ν1

ν2

ν3

ν4

ν5

ν6

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

2
2
2
3
3
3
4
4
4
5
5
5
6
6
6

3
3
3
2
2
2
2
2
2
2
2
2
2
2
2

4
5
6
4
5
6
3
5
6
3
4
6
3
4
5

5
4
4
5
4
4
5
3
3
4
3
3
4
3
3

6
6
5
6
6
5
6
6
5
6
6
4
5
5
4

1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

1
1
1
2
2
2
2
2
2
3
3
3
3
3
3

2
2
2
1
1
1
1
1
1
1
1
1
1
1
1

2
3
3
2
3
3
2
3
3
2
2
3
2
2
3

3
2
2
3
2
2
3
2
2
2
2
2
2
2
2

3
3
3
3
3
3
3
3
3
3
3
2
3
3
2

Table 2 The permutations ν ∈ d(P3 ) of the multi-set {1, 1, 2, 2, 3, 3}, their multiplicities #d −1 (ν) and cycle structures
No.

ν1

ν2

ν3

ν4

ν5

ν6

#d −1 (ν)

cycles

1
2
3
4
5
6

1
1
1
1
1
1

1
1
2
2
3
3

2
2
1
1
1
1

2
3
2
3
2
3

3
2
3
2
2
2

3
3
3
3
3
2

1
2
2
4
4
2

(1), (2), (3)
(1), (2,3)
(1,2), (3)
(1,2,3)
(1,3,2)
(1,3), (2)

multiplicities taken into account) as
E(ϕ(t)r ) =

∑
ν∈d(Pr )

#d −1 (ν)

Z

n

r

∑ ∏ K`ν2s−1 `ν2s (tν2s−1 − tν2s )dt1 × . . . × dtr .

[0,t]r ` ,...,`r =1 s=1
1

(126)

Now, any given permutation ν ∈ d(Pr ) can be partitioned into cycles c which pass through elements of the set
{1, . . . , r} twice. Every such cycle of period p can be represented as a p-tuple (c1 , c2 , . . . , c p ) of pairwise different

20

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

1)

1

1

2

2

3

3

2)

1

1

2

3

2

3

3)

1

2

1

2

3

3

4)

1

2

1

3

2

3

5)

1

3

1

2

2

3

6)

1

3

1

3

2

2

Fig. 1 These directed graphs represent the permutations ν of the multi-set {1, 1, 2, 2, 3, 3} in Tab. 2 (for the case r = 3) and their cycle structure.
The straight edges in each of the graphs depict the pairs (c j , c j+1 ) in (128) and are oriented rightwards (γ j = 0) or leftwards (γ j = 1), with the latter
indicating the inversions in these pairs. Each of the curved edges represents a transition to the remaining occurrence of the current element of the
cycle in the permutation. The fourth and fifth permutations are monocyclic. According to the second last column of Tab. 2, their total multiplicity is
8.

elements of the set {1, . . . , r} (see, for example, the rightmost column of Tab. 2). In terms of these cycles (which
depend on a particular permutation ν ∈ d(Pr )), the product in (126) takes the form
p

r

[γ ]

∏ K`ν2s−1 `ν2s (tν2s−1 − tν2s ) = ∏ ∏ K`c jj `c j+1 (tc j − tc j+1 ).

(127)

c j=1

s=1

Here, for a given cycle c := (c1 , . . . , c p ) of period p, the convention c p+1 := c1 is used, and γ1 , . . . , γ p are auxiliary
variables

0 if c j < c j+1
γ j := I(c j , c j+1 ) =
(128)
1 if c j > c j+1
which, in accordance with (114), indicate inversions for consecutive elements of the cycle c (we let γ1 := 0 in the case
p = 1); see Fig. 1. Also, use is made of the notation
(
K`c j `c j+1 (tc j − tc j+1 ) if c j < c j+1
[γ j ]
K`c `c (tc j − tc j+1 ) :=
.
(129)
K`c j+1 `c j (tc j+1 − tc j ) if c j+1 < c j
j j+1
Since the cycles of a given permutation ν partition the set {1, . . . , r}, then, in view of (128) and (129), the summation
of (127) over the independent indices `1 , . . . , `r leads to
n

∑

∏ K`ν2s−1 `ν2s (tν2s−1 − tν2s ) =

`1 ,...,`r =1 s=1

p

n

r

[γ ]

∑ ∏ ∏ K`c jj `c j+1 (tc j − tc j+1 )

`1 ,...,`r =1 c j=1
p
−
→ [γ ]
j

(tc j − tc j+1 ),

= ∏ Tr ∏ K
c

(130)

j=1

where the last equality also employs (111) and (112). By a similar reasoning, the integration of the right-hand side of
(130) over t1 , . . . ,tr yields
p

Z

Z
−
→ [γ ]
j (t − t
Tr
K
)dt
×
.
.
.
×
dt
=
cj
c j+1
r
1
∏ ∏
∏

[0,t]r c

j=1

c

[0,t] p

p

−
→
Tr ∏ K [γ j ] (t j − t j+1 )dt1 × . . . × dt p ,
j=1

(131)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

21

where t p+1 := t1 . Substitution of (131) into (126) leads to
p

r

E(ϕ(t) ) =

#d

∑

−1

−
→

Z

(ν) ∏
c

ν∈d(Pr )

[0,t] p

Tr ∏ K [γ j ] (t j − t j+1 )dt1 × . . . × dt p .

(132)

j=1

A combination of (132) with (62) and tedious combinatorial considerations (omitted here for brevity) show that the
contributions to the rth cumulant Kr (ϕ(t)) from all those permutations ν, which consist of at least two cycles, cancel
out in the sense that
p

#d −1 (ν) ∏

∑

Z

c

ν∈d(Pr ) with at least two cycles
r

= r! ∑

[0,t] p

−
→
Tr ∏ K [γ j ] (t j − t j+1 )dt1 × . . . × dt p
j=1

k

(−1)k
k

k=2

∑

j1 ,..., jk >1: j1 +...+ jk =r

E(ϕ(t) js )
.
js !
s=1

∏

Hence, only monocyclic permutations ν ∈ d(Pr ), consisting of a single cycle of period r (the set of such permutations
is denoted by Mr ), contribute to the cumulant:
r

Kr (ϕ(t)) =

∑

#d

−1

−
→

Z

(ν)

ν∈Mr

[0,t]r

Tr ∏ K [γ j ] (t j − t j+1 )dt1 × . . . × dtr .

(133)

j=1

The cycle, associated with a monocyclic permutation ν ∈ Mr , has the form (1, c2 , . . . , cr ), where (c2 , . . . , cr ) is an
arbitrary permutation of the set {2, . . . , r}. For any such (c2 , . . . , cr ), all monocyclic permutations ν ∈ Mr with the
cycle (1, c2 , . . . , cr ) are described by
ν = (1, c2 , 1, cr , p1 , . . . , pr−2 ),
(134)
where p1 , . . . , pr−2 are pairs of integers which are obtained by ordering the pairs
(a j , b j ) := (min(c j , c j+1 ), max(c j , c j+1 ))

(135)

for all j = 2, . . . , r − 1, so that their first entries a2 , . . . , ar−1 are ordered in ascending order according to (125). The
ordering is delivered by all those permutations f1 , . . . , fr−2 of the set {2, . . . , r − 1} which satisfy
a f1 6 . . . 6 a fr−2

(136)

and the trailing pairs in (134) are given in terms of (135) by
j = 1, . . . , r − 2.

p j := (a f j , b f j ),

(137)

Such a permutation f1 , . . . , fr−2 is unique only if a2 , . . . , ar−1 are pairwise different. If there are µ different pairs
of identical elements among a2 , . . . , ar−1 , then the number of permutations f1 , . . . , fr−2 which secure (136) is 2µ .
However, in both scenarios, there are 2r−1 permutations κ ∈ Pr such that the corresponding ν in (124) is monocyclic
→
−r
with the given cycle (1, c2 , . . . , cr ). Now, since the matrix ∏ j=1 K [γ j ] (t j −t j+1 ) on the right-hand side of (133) depends
on a monocyclic permutation ν in (134)–(137) only through the inversion indicators γ1 , . . . , γr in (128), which are
completely specified by c2 , . . . , cr , then
r

Kr (ϕ(t)) = 2r−1

Z

∑

r
(c2 ,...,cr ) [0,t]

−
→
Tr ∏ K [γ j ] (t j − t j+1 )dt1 × . . . × dtr ,

(138)

j=1

where the sum is over all (r −1)! permutations (c2 , . . . , cr ) of the set {2, . . . , r}. The summation in (138) can be reduced
to that over the binary variables γ2 , . . . , γr−1 as
r−1

Kr (ϕ(t)) = 2

Z

∑

γ2 ,...,γr−1 =0,1

∆r,γ2 ,...,γr−1

[0,t]r

r−1


−
→
Tr K(t1 − t2 ) ∏ K [γ j ] (t j − t j+1 )K(t1 − tr )T dt1 × . . . × dtr ,

(139)

j=2

since γ1 = 0 and γr = 1 in view of (128) and c1 = 1 < c2 and cr > cr+1 = 1. Here, in accordance with (115), the
coefficient ∆r,γ2 ,...,γr−1 is the number of those permutations (c2 , . . . , cr ) of the set {2, . . . , r} which have the inversion

22

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

indicators γ2 , . . . , γr−1 . Therefore, the sum of these coefficients is given by (116). From (122) and the symmetry of the
matrix Π , it follows that
r

r

j=1

j=1

√ 
−
→
−
→ √
Tr ∏ K [γ j ] (t j − t j+1 ) = Tr ∏ Π S[γ j ] (t j − t j+1 ) Π
r−1


−
→

= Tr Π S(t1 − t2 ) ∏ Π S[γ j ] (t j − t j+1 ) Π S(t1 − tr )T .

(140)

j=2

Substitution of (140) into (139) leads to (117), thus completing the proof.

Although Theorem 3 is formulated and proved for the case Π < 0, the relation (117) remains valid irrespective of
the matrix Π being positive semi-definite.
In contrast to covariance functions of classical random processes, the quantum covariance function S in (102) is
not invariant under the transformation (111) since S[1] (τ) = S(τ) 6= S(τ) in general. Therefore, the integral in (117)
has a more complicated structure than the r-fold convolution of the function Π S with itself. For example, in the case
r = 3, the sum in (117) is over a binary index and consists of two terms. In view of (116) and the invariance of the
matrix trace with respect to the transpose and cyclic permutations, this leads to
Z


K3 (ϕ(t)) =4 ∆3,0
Tr Π S(t1 − t2 )Π S(t2 − t3 )Π S(t1 − t3 )T dt1 dt2 dt3
[0,t]3

Z

+ ∆3,1

[0,t]3



Tr Π S(t1 − t2 )Π S(t3 − t2 )T Π S(t1 − t3 )T dt1 dt2 dt3
Z

=4(∆3,0 + ∆3,1 )
Z

=8

[0,t]3

[0,t]3


Tr Π S(t1 − t2 )Π S(t2 − t3 )Π S(t1 − t3 )T dt1 dt2 dt3

Π , S(t1 − t2 )Π S(t2 − t3 )Π S(t1 − t3 )T dt1 dt2 dt3 .

(141)

For larger values of r, the dependence of the coefficient ∆r,γ in (117) on the multi-index γ ∈ {0, 1}r−2 is fairly complicated and manifests fractal-like fluctuations, as illustrated by Fig. 2.
The following corollary of Theorem 3 describes the infinite-horizon asymptotic behaviour of the cumulants of
arbitrary order, including its frequency-domain representation. To this end, we denote by D the spectral density of the
quantum covariance function (102):
Z +∞

D(λ ) :=

S(τ)e−iλ τ dτ.

(142)

−∞

Due to the Bochner-Khinchin criterion of positive semi-definiteness of Hermitian kernels [16] and the standard properties of the Fourier transforms, D(λ ) = D(λ )∗ < 0 for all frequencies λ ∈ R. Furthermore, D admits the factorization
D(λ ) = G(iλ )Ω G(iλ )∗

(143)

in terms of the quantum Ito matrix Ω from (10) and the transfer function G from W to X in (9) associated with the
matrix pair (A, B) by
Z +∞

G(s) :=
0

e−st e−tA Bdt = (sIn − A)−1 B

(144)

in the closed right half-plane Res > 0 (since A is Hurwitz). The transformation (111), applied to the covariance function
S, carries over to the spectral density D in (142) as
S[1] (τ) = S(−τ)T
Z
1 +∞ −iλ τ
e
D(λ )T dλ
=
2π −∞
Z
Z
1 +∞ iλ τ
1 +∞ iλ τ [1]
T
=
e D(−λ ) dλ =
e D (λ )dλ ,
2π −∞
2π −∞
with the latter being the inverse Fourier transform of D[1] . Also, similarly to (112), for appropriately dimensioned
matrix-valued functions f1 , . . . , fr on the real line, their convolution is denoted by
r

→
−
~ fk := f1 ∗ . . . ∗ fr .

k=1

(145)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

23

#10 4

5

4.5

4

"r, .

2

,..., .

r-1

3.5

3

2.5

2

1.5

1

0.5

50

100

150

200

250

'r-3
. 2j,
j=0 j+2

300

350

400

450

500

r=11

Fig. 2 The dependence of the coefficient ∆r,γ in (117) on the binary (r − 2)-index γ in the case r = 11. The multi-index encodes the integers
0, . . . , 2r−2 − 1 over the horizontal axis of the graph.

Theorem 4 Under the conditions of Theorem 3, for any r > 2, the asymptotic growth rate of the rth cumulant in (117)
can be computed as

lim

t→+∞

1
t


Kr (ϕ(t)) = 2r−1

∑

r−1



→
−
∆r,γ Tr Π S ∗ ~ Π S[γ j ] ∗ Π S[1] (0)
j=2

γ∈{0,1}r−2

2r−2
=
∑ ∆r,γ
π γ∈{0,1}
r−2

Z +∞
−∞

r−1


−
→

Tr Π D(λ ) ∏ Π D[γ j ] (λ ) Π D[1] (λ ) dλ ,

(146)

j=2

where the notation (145) is used, and D is the spectral density of the steady-state quantum covariance function S in
(102) given by (142)–(144).

Proof The relations (146) are obtained by applying Lemma 6 of Appendix C to the integrals in (117) with the functions
f1 := Π S, fr := Π S[1] and fk := Π S[γk ] for k = 2, . . . , r − 1 whose Fourier transforms are Π D, Π D[1] and Π D[γk ] ,
respectively.

Note that in the case r = 2, the first of the equalities (146) reproduces (106). In application to the third cumulant
in (141), Theorem 4 yields
lim

t→+∞

1
t



K3 (ϕ(t)) = 8Tr Π S ∗ Π S ∗ Π S[1] (0)
Z

=8
=

4
π

Π , S(σ )Π S(τ − σ )Π S(τ)T dσ dτ

R2
Z +∞
−∞

Π , (D(λ )Π )2 D(−λ )T dλ .

24

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

In comparison with the classical case, the representation (146) involves the sum of terms whose number grows exponentially with r. The complicated structure of these terms, which are controlled by successive inversions in permutations, and the combinatorial nature of the coefficients ∆r,γ give rise to the problem of developing a recurrence relation
for the cumulant rates.

7 Large deviations estimates
Theorems 3 and 4 of the previous section can be employed for obtaining large deviations estimates for the process ϕ
in (50).3 More precisely, for any given time t > 0, application of the Cramer inequality [57] (see also [12, Section 3.5])
to the probability distribution Et of ϕ(t) yields the inequality
1

1
ln Et ([εt, +∞)) 6 inf
ln Ξθ (t) − θ ε
θ >0 t
t

(147)

in terms of the QEF (52) for any ε > 0. Similarly to the classical case, the function being minimised in (147) is convex
with respect to θ and vanishes at θ = 0. The following theorem leads to an upper bound for the right-hand side of
(147)4 , uniform in time t. For its formulation, we need an auxiliary function N associated with the quantum covariance
function S in (102) by
√
√
(148)
N(τ) := k Π S(τ) Π k,
where k · k is the `2 -induced operator norm of a matrix. Note that N is an R+ -valued integrable function which is
even in view of the second equality in (45) and the invariance of the matrix norm k · k under the complex conjugate
transpose:
√
√
N(−τ) = k( Π S(τ) Π )∗ k = N(τ).
(149)
Theorem 5 Suppose the conditions of Theorem 3 are satisfied. Also, let the risk-sensitivity parameter θ belong to the
interval
1
06θ <
,
(150)
2kFk∞
where kFk∞ is the L∞ -norm of the Fourier transform F of the function N from (148) given by F(λ ) :=
Then the QEF (52) admits the upper bound
sup

R +∞
−∞

N(τ)e−iλ τ dλ .

Z

n +∞
ln Ξθ (t) 6 −
ln(1 − 2θ F(λ ))dλ .
t
4π −∞

1

t>0

(151)


Proof Similarly to (149), from the invariance of the operator norm k · k with respect to the transpose and the symmetry
of the function N, it follows that
√
√
√
√
k Π S[1] (τ) Π k = k( Π S(−τ) Π )T k = N(τ),
(152)
where use is made of the transformation (111). A combination of (152) with the submultiplicativity of the operator
norm leads to
Z
[0,t]r

r−1


−
→

Tr Π S(t1 − t2 ) ∏ Π S[γ j ] (t j − t j+1 ) Π S(t1 − tr )T dt1 × . . . × dtr
j=2

r−1

Z

6n

[0,t]r

Z

6n

√
√
√ −
√ √
→ √
Π S(t1 − t2 ) Π ∏ Π S[γ j ] (t j − t j+1 ) Π Π S(t1 − tr )T Π dt1 × . . . × dtr
j=2

r

∏ N(t j − t j+1 )dt1 × . . . × dtr

[0,t]r j=1

6 ntN ∗r (0) =
3
4

nt
2π

Z +∞

F(λ )r dλ ,

(153)

−∞

Such estimates are physically meaningful in the case Π < 0 when the self-adjoint quantum variable ϕ(t) is a positive semi-definite operator.
which is the negative of the Legendre transform of 1t ln Ξθ (t) as a function of θ > 0

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

25

where N ∗r denotes the r-fold convolution of the function N with itself. The last inequality in (153) follows from the
nonnegativeness of N and the identity (C7) in the proof of Lemma 6, while the first equality uses the upper bound
|TrK| 6 nkKk which holds for any matrix K ∈ Cn×n . In view of (153) and (116), the rth cumulant in (117) satisfies
the inequality
Kr (ϕ(t)) 6 2r−1

∑

∆r,γ

γ∈{0,1}r−2

= (r − 1)!2r−2

nt
π

nt
2π

Z +∞

Z +∞

F(λ )r dλ

−∞

F(λ )r dλ

(154)

−∞

for any r = 1, 2, 3, . . .. Assuming that θ is nonnegative and sufficiently small in the sense of (150), substitution of (154)
into (60) leads to
+∞

ln Ξθ (t) 6
=

nt
θr
∑ r! (r − 1)!2r−2 π
r=1
nt
4π

+∞

1
∑r
r=1

nt
=−
4π

Z +∞

Z +∞

F(λ )r dλ

−∞

(2θ F(λ ))r dλ

−∞

Z +∞

ln(1 − 2θ F(λ ))dλ .

(155)

−∞

The inequality (151) can now be obtained by dividing both sides of (155) by t > 0 and using the fact that t is otherwise
arbitrary.

A combination of (147) with (151) of Theorem 5 leads to
sup
t>0

1
t


ln Et ([εt, +∞)) 6



inf

1
06θ < 2kFk

−

n
4π

Z +∞
−∞


ln(1 − 2θ F(λ ))dλ − θ ε .

(156)

∞

The function under minimization in (156) is strictly convex with respect to θ in the interval (150) and hence, has a
unique minimum. The optimal value of θ can be found from the equation
−

n
∂θ
4π

Z +∞

ln(1 − 2θ F(λ ))dλ =

−∞

n
2π

Z +∞
−∞

F(λ )
dλ = ε
1 − 2θ F(λ )

(157)

whose left-hand side is a strictly increasing function of θ . The solution θ > 0 exists and is unique if the “scale”
parameter ε is large enough in the sense that
ε
1
>
n 2π

Z +∞
−∞

√
√
F(λ )dλ = N(0) = k Π (P + iΘ ) Π k.

(158)

Note that (151), (156) and (157) remain valid if the function N in (148) is replaced with its even upper estimate. Since
the quantum covariance function S decays exponentially fast at infinity, such an upper bound for N can be found among
functions with a rational Fourier transform F. In this case, the integral in (156) lends itself to effective evaluation [44]
as does the solution of (157).
Theorem 6 Suppose the conditions of Theorem 3 are satisfied. Then the tail of the probability distribution Et of the
process ϕ(t) in (50) admits an upper bound
sup
t>0

where

1

 nµ 
nα
ε 
ln Et ([εt, +∞)) 6
2−
−
,
t
4
ε
nα
√ √
√
α := k Π Γ kkΓ −1/2 (P + iΘ ) Π k.

ε > nα,

(159)

(160)

Here, (µ,Γ ) is any pair of a scalar µ > 0 and a real positive definite symmetric matrix Γ of order n satisfying the
algebraic Lyapunov inequality (ALI)
AΓ + Γ AT 4 −2µΓ .

(161)

26

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

Proof The pairs (µ,Γ ), satisfying (161), exist due to the matrix A being Hurwitz. Any such pair satisfies the differential matrix inequality
T
T
∂τ e2µτ eτAΓ eτA = e2µτ eτA (AΓ + Γ AT + 2µΓ )eτA 4 0.
(162)
T

As a matrix-valued version of the Gronwall-Bellman lemma, the integration of (162) over τ leads to eτAΓ eτA 4
e−2µτ Γ , which is equivalent to the contraction property
√
kΓ −1/2 eτA Γ k 6 e−µτ ,
τ > 0.
(163)
A combination of (163) with (102) and the submultiplicativity of the operator norm allows the exponential decay of
the function (148) to be quantified as
√ √
√
√
N(τ) = k Π Γ Γ −1/2 eτA Γ Γ −1/2 (P + iΘ ) Π k
√
√
√ √
6 k Π Γ kkΓ −1/2 eτA Γ kkΓ −1/2 (P + iΘ ) Π k
6 αe−µτ

(164)

for all τ > 0, where α is given by (160). Since the function N is even, the inequality (164) can be extended to all τ as
b
N(τ) 6 αe−µ|τ| =: N(τ),

τ ∈ R.

(165)

The right-hand side of (165) has the Fourier transform
Z +∞

b ) :=
F(λ

−∞

−iλ τ
b
dτ =
N(τ)e

2α µ
λ 2 + µ2

,

(166)

which is a rational function of λ with the L∞ -norm
b ∞=
kFk

2α
.
µ

(167)

The corresponding integral in (156) can be evaluated by using a contour integration or through the parametric differentiation
−∂θ

Z +∞
−∞

b ))dλ = 2
ln(1 − 2θ F(λ

Z +∞
−∞

b )
F(λ
dλ
b )
1 − 2θ F(λ

Z +∞

= 4α µ
=p

−∞

1
dλ
λ 2 + µ 2 − 4θ α µ

4πα µ
µ 2 − 4θ α µ

.

(168)

After integration over θ , this yields
−

Z +∞
−∞

b ))dλ = 2π µ −
ln(1 − 2θ F(λ

p

µ 2 − 4θ α µ



(169)

for any θ satisfying
06θ <

1
µ
=
b
4α
2kFk∞

(170)

b and Fb in (165) and (166), and
in view of (150) and (167). By appropriately modifying (156)–(158) for the functions N
using (168) and (169), it follows that

 n Z +∞

1
b ))dλ − θ ε
sup
ln Et ([εt, +∞)) 6
inf
−
ln(1 − 2θ F(λ
4π −∞
06θ < b1
t>0 t
2kFk∞

n
p

= inf µ
µ − µ 2 − 4θ α µ − θ ε
06θ < 4α 2
 nα 2 
nα  µε 
nµ 
=
1−
−
1−
2
ε
4α
ε
nµ 
nα
ε 
=
2−
−
,
4
ε
nα

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

27

µ
2
where the minimum over the interval (170) is achieved at θ = 4α
(1 − ( nα
ε ) ), provided ε > nα, thus establishing
(159).

Note that the ALI (161) has a positive definite solution Γ for any µ < − max{Reλ : λ ∈ S}, where S is the
spectrum of the Hurwitz matrix A. A particular choice of such a matrix Γ (which can be carried out, for example,
by using the eigenbasis of A in the case when A is diagonalizable) affects the constant α in (160) which enters (159)
together with µ.

8 Classical covariance correspondence
We will now discuss a correspondence between vectors of self-adjoint quantum variables and auxiliary classical random vectors at the level of the second-order moments. More precisely, let ζ be a classical Cn -valued random vector
given by
ζ := ξ + iη,
ξ := Reζ ,
η := Imζ .
(171)
Here, ξ and η are classical random vectors with values in Rn , which are also assembled into the R2n -valued vector
 
ξ
ϑ :=
.
(172)
η
Suppose ζ has zero mean and is square integrable, that is, Eζ = 0 and E(|ζ |2 ) < +∞, where |ζ |2 = |ξ |2 + |η|2 . Then
its covariance matrix is given by
 

 
1
cov(ζ ) = E(ζ ζ ∗ ) = ( 1 i ⊗ In )cov(ϑ )
⊗ In
−i
= E(ξ ξ T ) + E(ηη T ) − i(E(ξ η T ) − E(ηξ T ))
= cov(ξ ) + cov(η) − 2iA(cov(ξ , η)),

(173)

where A(M) := 21 (M − M T ) denotes the antisymmetrizer of square matrices which is applied here to the crosscovariance matrix cov(ξ , η) = E(ξ η T ) of the zero-mean vectors ξ and η. Of particular interest for our purposes
is the case when ξ and η in (172) have identical covariance matrices and an antisymmetric cross-covariance matrix:




1 P −Θ
E(ξ ξ T ) E(ξ η T )
cov(ϑ ) =
=
.
(174)
E(ηξ T ) E(ηη T )
2 Θ P
Here, P is a real positive semi-definite symmetric matrix of order n, and Θ is a real antisymmetric matrix of order n.
In this case, (173) takes the form
cov(ζ ) = P + iΘ ,
(175)
h
i
P −Θ
where P + iΘ is isospectral (up to multiplicity of eigenvalues) to the matrix Θ P on the right-hand side of (174),
with both matrices being positive semi-definite. Now, let ζ depend on time and form a Cn -valued Markov diffusion
process governed by a classical SDE
1
dζ = Aζ dt + √ BΩ dω
(176)
2
driven by a standard Wiener process ω in Rm . Here, Ω is the matrix from (10), and the matrices A and B are given by
(13) as before. Since both A and B are real, then, due to the structure of Ω , the SDE (176) can be represented in terms
of the real and imaginary parts in (171) as
1
dξ = Aξ dt + √ Bdω,
2

1
dη = Aηdt + √ BJdω.
2

(177)

Equivalently, the augmented process ϑ in (172) satisfies the SDE
 
1
I
dϑ = (I2 ⊗ A)ϑ dt + √ (I2 ⊗ B) m dω,
J
2

(178)

The following lemma establishes a correspondence between the classical Markov diffusion processes ξ , η, ζ , and the
linear quantum stochastic system considered in the previous sections.

28

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

Lemma 3 Suppose the matrices A and B are given by (13), with A being Hurwitz. Then the classical Markov diffusion
process ϑ in (172), governed by (178), has a unique invariant measure (the invariant joint distribution of ξ and η
in (177)) which is a Gaussian distribution in R2n with zero mean and covariance matrix (174), where P is given by
(24) and Θ is the CCR matrix from (5). The corresponding invariant measure of the process ζ in (171) and (176) is
a Gaussian distribution in Cn with zero mean and the covariance matrix (175). Furthermore, if ϑ and the quantum
process X are initialised at the corresponding invariant Gaussian states, then their two-point covariance functions
reproduce each other in the sense that
E(ζ (t)ζ (s)∗ ) = E(X(t)X(s)T ) = e(t−s)A (P + iΘ )

(179)

for all t > s > 0, with the classical expectation on the left-hand side and the quantum expectation on the right-hand
side.

Proof Due to A being Hurwitz, the existence, uniqueness and Gaussian nature of the invariant measure for the process
ϑ follows from (178), whereby the invariant distribution has zero mean and covariance matrix cov(ϑ ) satisfying the
ALE


1
Im −J
T
(I2 ⊗ BT ) = 0,
(180)
(I2 ⊗ A)cov(ϑ ) + cov(ϑ )(I2 ⊗ A ) + (I2 ⊗ B)
J Im
2
where use is made of the orthogonality of the real antisymmetric matrix J in (11). The ALE (180) splits into three
equations
1
Acov(ξ ) + cov(ξ )AT + BBT = 0,
2
1
Acov(η) + cov(η)AT + BBT = 0,
2
1
Acov(ξ , η) + cov(ξ , η)AT − BJBT = 0.
2

(181)
(182)
(183)

Since the ALEs (181) and (182) are identical to (25) up to a factor of 12 (and have unique solutions due to the matrix
A being Hurwitz), it follows that
1
cov(ξ ) = cov(η) = P,
(184)
2
with P given by (24). By a similar reasoning, a comparison of (183) with the PR condition (14) leads to
1
cov(ξ , η) = − Θ ,
2

(185)

where Θ is the CCR matrix from (5). The relations (184) and (185) imply that the covariance matrix of the invariant
joint Gaussian distribution of the processes ξ and η is indeed given by (174), with (175) describing the corresponding
covariance matrix for ζ . Now, the linear SDE (176) implies that
1
ζ (t) = e(t−s)A ζ (s) + √
2

Z t

e(t−τ)A BΩ dω(τ)

(186)

s

for all t > s > 0. Due to ω being a standard Wiener process independent of the initial value ζ (0), it follows from (186)
that
E(ζ (t)ζ (s)∗ ) = e(t−s)A E(ζ (s)ζ (s)∗ ).
(187)
If the quantum system is initialised at the invariant Gaussian state and ϑ (0) has the corresponding invariant Gaussian
distribution, then
cov(ζ (s)) = P + iΘ = cov(X(s))
(188)
for all s > 0. Substitution of (188) into (48) and (187) leads to (179).

Note that Lemma 3 is concerned only with the correspondence between classical Gaussian Markov diffusion
processes and quantum system variables in Gaussian quantum states at the level of covariances. This correspondence
involves complex conjugation and does not extend, in general, from the second-order moments to arbitrary higherorder moments. For example,
ζ ∗Π ζ = ξ TΠ ξ + η TΠ η
(189)
is a real-valued classical random process which has the same steady-state mean value as the quantum process ψ in
(99):
E(ζ ∗ Π ζ ) = hΠ , E(ζ ζ ∗ )i = hΠ , P + iΘ i = hΠ , Pi,

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

29

which follows from (179). At the same time, application of the classical covariance relation (108) to (189) leads to
var(ζ ∗ Π ζ ) = var(ξ T Π ξ ) + var(η T Π η) + 2cov(ξ T Π ξ , η T Π η)
= 2 hΠ , cov(ξ )Π cov(ξ ) + cov(η)Π cov(η) + 2cov(ξ , η)Π cov(η, ξ )i
= hΠ , PΠ P −Θ ΠΘ i ,

(190)

where use is also made of (184) and (185). Note that the right-hand side of (190) is different from its quantum
counterpart
var(X T Π X) = 2 hΠ , PΠ P +Θ ΠΘ i
which is obtained by letting σ = τ in (101) and (102). This discrepancy (which manifests itself already at the level
of the fourth-order moments) is closely related to the absence of classical joint probability distributions for noncommuting quantum variables. Nevertheless, of interest is further comparison of the QEF (52) with the corresponding
risk-sensitive cost functional for the SDE (176) whose asymptotic behaviour is described by
Z
1

Rt
∗
1 +∞
ln det(In − θ Π G(iλ )Ω G(iλ )∗ )dλ
ln Eeθ 0 ζ (s) Π ζ (s)ds = −
lim
t→+∞ t
4π −∞
Z

1 +∞ θ r +∞
=
Tr (Π D(λ ))r dλ .
(191)
∑
4π r=1 r −∞
Here, G is the transfer function in (144), so that 21 G(iλ )Ω G(iλ )∗ is the spectral density of the process ζ in (176)
which is expressed in terms of the spectral density D for the quantum covariance function S in (142) and (143). In the
case Π < 0, the right-hand side of (191) is well-defined if the risk sensitivity parameter θ > 0 is bounded in terms of
2
, where we have also used the property Ω 2 = 2Ω
a weighted H∞ -norm of the transfer function G as θ < k√Π GΩ
k2∞

R +∞
of the matrix Ω from (10). Note that −∞ Tr (Π D(λ ))r dλ in (191) is a classical counterpart of the weighted sum of
integrals on the right-hand side of (146).
9 A numerical example of the quartic approximation of the QEF and large deviations bounds
We will now demonstrate the computation of the quartic approximation of the QEF from Section 5 (in the infinitehorizon limit) for a two-mode OQHO with n = 4 driven by m = 4 external fields. It is assumed that the CCR matrix in
(5) is given by
1
Θ = J ⊗ I2 ,
(192)
2
which, in accordance with (6) and (7), corresponds to the system variables q1 , q2 , p1 , p2 consisting of the conjugate
quantum mechanical positions and momenta. The energy and coupling matrices of the OQHO were randomly generated so as to make the matrix A in (13) Hurwitz:




−0.1027 1.3449 −0.2403 −1.3994

1.3449 1.5008 −0.1856 0.9212
R = −0.2403 −0.1856 −0.5704 −0.4146,
−1.3994 0.9212 −0.4146 −0.3233

0.8726

0.1632

2.1844 −1.9270

0.1179 −0.8147 −0.0938 0.5214 
M=
−1.5031 0.4037 −0.2942 −2.0544
0.9218 0.7562 −0.5048 −0.2698

(193)

(the resulting A has the eigenvalues −0.5532 ± 2.5929i, −1.3302, −4.2068). The positive definite weighting matrix
Π , specifying the QEF through (50)–(52), was also randomly generated:


3.5050 −0.5447 0.0672 −2.3918

−0.5447 4.0758 −1.1876 0.0215
Π =  0.0672 −1.1876 5.1422 −1.4628.

(194)

−2.3918 0.0215 −1.4628 4.5416

The infinite-horizon controllability Gramian P in (24) and the matrix T in (107) were found by successively solving
the ALEs (25) and (97):




3.7981 −2.5143 −3.8716 −1.6214

−2.5143 4.9443 0.5356 0.4305
P = −3.8716 0.5356 6.7086 2.8509 ,
−1.6214 0.4305

2.8509

1.4473

131.5431 −108.9564 −138.4442 −58.4033

−108.9564 138.7545 60.4808 21.2105
T = −138.4442 60.4808 204.6153 91.4998 .
−58.4033

21.2105

91.4998

41.2158

The corresponding asymptotic growth rates for the first two cumulants of the quantum process ϕ in (94) and (96) took
the values

1

1
lim
Eϕ(t) = 74.9147,
lim
var(ϕ(t)) = 8.9399 × 103 .
t→+∞ t
t→+∞ t

30

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

In this example, the threshold value for the risk-sensitivity parameter θ in (110), within which the quartic term in (93)
can be neglected, is
θ0 = 0.0168.
k

Since θ0 is small, the coefficients θk! of the higher-order moments E(ϕ(t)k ) and cumulants Kk (ϕ(t)) in (59) and
(60) (with k > 3) are negligible compared to θ 2 for such values of θ . This makes the quartic approximation (93) a
reasonable approximation of the QEF for the range of small values of θ in this example. Since the matrix A under
consideration is diagonalizable, the ALI (161) is satisfied with µ = 0.5532 (the negative of the largest real part of its
eigenvalues) and


1.4750 −0.4852 −1.4090 −0.2636

−0.4852 0.6271 0.2354 0.1475
Γ = UU = −1.4090 0.2354 1.6303 0.3569 ,
∗

−0.2636 0.1475

0.3569

0.2676

where the matrix U is formed from the eigenvectors of A. The corresponding constant (160) is α = 69.6784, and the
large deviations bound (159) of Theorem 6 is shown in Fig. 3.

0

n7(2-n,/0-0/(n,))/4

-0.5

-1

-1.5

-2

-2.5

400

600

800

1000

1200

1400

1600

1800

0

Fig. 3 The large deviations probability bound (159) as a function of the scale parameter ε > nα for the two-mode OQHO of the example specified
by (192)–(194).

10 Conclusion
For an open quantum harmonic oscillator, governed by a linear QSDE driven by vacuum bosonic fields, we have
considered a quadratic-exponential functional which penalizes the second and higher-order moments of the system
variables. We have obtained an integro-differential equation for the time evolution of the QEF and compared it with
the original quantum risk-sensitive performance criterion which was used previously in measurement-based quantum
control and filtering problems. We have discussed multi-point Gaussian quantum states for the system variables at different instants and employed their first four moments for a quartic approximation of the QEF whose infinite-horizon
asymptotic behaviour has also been investigated. A numerical example has been provided in order to demonstrate the
approximation. Higher-order cumulants, associated with the QEF, have been related to combinatorics of consecutive
inversions in permutations and gives rise to a nontrivial problem of their recursive computation. They have also been
used for a Cramer-type large deviations estimate for the system variables. We have also considered an auxiliary classical Gaussian Markov diffusion process in a complex Euclidean space, which reproduces the quantum system variables

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

31

at the level of covariance functions but has different higher-order moments relevant to the risk-sensitive criteria. The
results of the paper may find applications to coherent (measurement-free) quantum risk-sensitive control problems,
where the plant and controller form a fully quantum closed-loop system, and other settings with nonquadratic cost
functionals.

References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.

S.L.Adler, Derivation of the Lindblad operator structure by use of the Ito stochastic calculus, Phys. Lett. A, vol. 26, 2000, pp. 58–61.
B.D.O.Anderson, and J.B.Moore, Optimal Filtering, Prentice Hall, New York, 1979.
V.P.Belavkin, On the theory of controlling observable quantum systems, Autom. Rem. Contr., vol. 44, no. 2, 1983, pp. 178–188.
V.P.Belavkin, Noncommutative dynamics and generalized master equations, Math. Notes, vol. 87, no. 5, 2010, pp. 636–653.
A.Bensoussan, and J.H.van Schuppen, Optimal control of partially observable stochastic systems with an exponential-of-integral performance
index, SIAM J. Control Optim., vol. 23, 1985, pp. 599–613.
P.Billingsley, Convergence of Probability Measures, John Wiley & Sons, New York, 1968.
L.Bouten, and R.van Handel, On the separation principle of quantum control, arXiv:math-ph/0511021v2, August 22, 2006.
L.Bouten, R.Van Handel, M.R.James, An introduction to quantum filtering, SIAM J. Control Optim., vol. 46, no. 6, 2007, pp. 2199–2241.
D.J.Cross, and R.Gilmore, A Schwinger disentangling theorem, J. Math. Phys., vol. 51, 2010, pp. 103515-1–5.
C.D.Cushen, and R.L.Hudson, A quantum-mechanical central limit theorem, J. Appl. Prob., vol. 8, no. 3, 1971, pp. 454–469.
C.D’Helon, A.C.Doherty, M.R.James, and S.D.Wilson, Quantum risk-sensitive control, Proc. 45th IEEE CDC, San Diego, CA, USA, December
13–15, 2006, pp. 3132–3137.
P.Dupuis, and R.S.Ellis, A Weak Convergence Approach to the Theory of Large Deviations, Wiley, New York, 1997.
S.C.Edwards, and V.P.Belavkin, Optimal quantum filtering and quantum feedback control, arXiv:quant-ph/0506018v2, August 1, 2005.
G.B.Folland, Harmonic Analysis in Phase Space, Princeton University Press, Princeton, 1989.
C.W.Gardiner, and P.Zoller, Quantum Noise. Springer, Berlin, 2004.
I.I.Gikhman, and A.V.Skorokhod, The Theory of Stochastic Processes, Springer, Berlin, 2004.
V.Gorini, A.Kossakowski, E.C.G.Sudarshan, Completely positive dynamical semigroups of N-level systems, J. Math. Phys., vol. 17, no. 5,
1976, pp. 821–825.
J.Gough, and M.R.James, Quantum feedback networks: Hamiltonian formulation, Commun. Math. Phys., vol. 287, 2009, pp. 1109–1132.
A.S.Holevo, Quantum stochastic calculus, J. Sov. Math., vol. 56, no. 5, 1991, pp. 2609–2624.
A.S.Holevo, Exponential formulae in quantum stochastic calculus, Proc. Roy. Soc. Edinburgh, vol. 126A, 1996, pp. 375–389.
A.S.Holevo, Statistical Structure of Quantum Theory, Springer, Berlin, 2001.
R.A.Horn, and C.R.Johnson, Matrix Analysis, Cambridge University Press, New York, 2007.
R.L.Hudson, and K.R.Parthasarathy, Quantum Ito’s formula and stochastic evolutions, Commun. Math. Phys., vol. 93, 1984, pp. 301–323
L.Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables,
Biometrika, vol. 12, 1918, pp. 134–139.
K.Jacobs, and P.L.Knight, Linear quantum trajectories: applications to continuous projection measurements, Phys. Rev. A, vol. 57, no. 4, 1998,
pp. 2301–2310.
D.H.Jacobson, Optimal stochastic linear systems with exponential performance criteria and their relation to deterministic differential games,
IEEE Trans. Autom. Control, vol. 18, 1973, pp. 124–31.
M.R.James, Risk-sensitive optimal control of quantum systems, Phys. Rev. A, vol. 69, 2004, pp. 032108-1–14.
M.R.James, A quantum Langevin formulation of risk-sensitive optimal control, J. Opt. B, vol. 7, 2005, pp. S198–S207.
M.R.James, and J.E.Gough, Quantum dissipative systems and feedback control design by interconnection, IEEE Trans. Automat. Contr., vol.
55, no. 8, 2008, pp. 1806–1821.
M.R.James, H.I.Nurdin, and I.R.Petersen, H ∞ control of linear quantum stochastic systems, IEEE Trans. Automat. Contr., vol. 53, no. 8, 2008,
pp. 1787–1803.
S.Janson, Gaussian Hilbert Spaces, Cambridge University Press, Cambridge, 1997.
I.Karatzas, and S.E.Shreve, Brownian Motion and Stochastic Calculus, 2nd Ed., Springer, New York, 1991.
T.Kato, Perturbation Theory for Linear Operators, 2nd Ed., Springer, Berlin, 1976.
H.Kwakernaak, and R.Sivan, Linear Optimal Control Systems, Wiley, New York, 1972.
L.D.Landau, and E.M.Lifshitz, Quantum Mechanics: Non-relativistic Theory, 3rd Ed., Pergamon Press, Oxford, 1991.
G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys., vol. 48, 1976, pp. 119–130.
M.Lutzky, Parameter differentiation of exponential operators and the Baker-Campbell-Hausdorff formula, J. Math. Phys., vol. 9, no. 7, 1968,
pp. 1125–1128.
A.I.Maalouf, and I.R.Petersen, Coherent LQG control for a class of linear complex quantum systems, IEEE European Control Conference,
Budapest, Hungary, 23-26 August 2009, pp. 2271–2276.
M.Maggiore, A Modern Introduction to Quantum Field Theory, Oxford University Press, New York, 2005.
W.Magnus, On the exponential solution of differential equations for a linear operator, Comm. Pure Appl. Math., vol. 7, no. 4, 1954, pp. 649–673.
J.R.Magnus, The moments of products of quadratic forms in normal variables, Statist. Neerland., vol. 32, 1978, pp. 201–210.
E.Merzbacher, Quantum Mechanics, 3rd Ed., Wiley, New York, 1998.
P.-A.Meyer, Quantum Probability for Probabilists, Springer, Berlin, 1995.
D.Mustafa, and K.Glover, Minimum Entropy H ∞ -Control, Springer Verlag, Berlin, 1990.
M.A.Nielsen, and I.L.Chuang, Quantum Computation and Quantum Information, Cambridge University Press, Cambridge, 2000.
H.I.Nurdin, M.R.James, and I.R.Petersen, Coherent quantum LQG control, Automatica, vol. 45, 2009, pp. 1837–1846.
K.R.Parthasarathy, An Introduction to Quantum Stochastic Calculus, Birkhäuser, Basel, 1992.
K.R.Parthasarathy, What is a Gaussian state? Commun. Stoch. Anal., vol. 4, no. 2, 2010, pp. 143–160.
K.R.Parthasarathy, and K.Schmidt, Positive Definite Kernels, Continuous Tensor Products, and Central Limit Theorems of Probability Theory,
Springer-Verlag, Berlin, 1972.
K.R.Parthasarathy, and R.Sengupta, From particle counting to Gaussian tomography, Inf. Dim. Anal., Quant. Prob. Rel. Topics, vol. 18, no. 4,
2015, pp. 1550023.

32

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

51. K.R.Parthasarathy, Quantum stochastic calculus and quantum Gaussian processes, Indian Journal of Pure and Applied Mathematics, 2015, vol.
46, no. 6, 2015, pp. 781–807.
52. I.R.Petersen, Guaranteed non-quadratic performance for quantum systems with nonlinear uncertainties, American Control Conference (ACC),
4-6 June 2014, arXiv:1402.2086 [quant-ph], 10 February 2014.
53. I.R.Petersen, Quantum linear systems theory, Open Automat. Contr. Syst. J., vol. 8, 2017, pp. 67–93.
54. I.R.Petersen, V.A.Ugrinovskii, and M.R.James, Robust stability of uncertain linear quantum systems, Phil. Trans. R. Soc. A, vol. 370, 2012,
pp. 5354–5363.
55. J.J.Sakurai, Modern Quantum Mechanics, Addison-Wesley, Reading, Mass., 1994.
56. A.J.Shaiju, and I.R.Petersen, A frequency domain condition for the physical realizability of linear quantum systems, IEEE Trans. Automat.
Contr., vol. 57, no. 8, 2012, pp. 2033–2044.
57. A.N.Shiryaev, Probability, 2nd Ed., Springer, New York, 1996.
58. V.S.Vladimirov. Methods of the Theory of Generalized Functions, London: Taylor & Francis, 2002.
59. I.G.Vladimirov, and I.R.Petersen, Hardy-Schatten norms of systems, output energy cumulants and linear quadro-quartic Gaussian control, Proc.
19th Int. Symp. Math. Theor. Networks Syst., Budapest, Hungary, July 5–9, 2010, pp. 2383–2390.
60. I.G.Vladimirov, and I.R.Petersen, A dynamic programming approach to finite-horizon coherent quantum LQG control, Australian Control
Conference, Melbourne, 10–11 November, 2011, pp. 357–362, arXiv:1105.1574v1 [quant-ph], 9 May 2011.
61. I.G.Vladimirov, and I.R.Petersen, A quasi-separation principle and Newton-like scheme for coherent quantum LQG control, Syst. Contr. Lett.,
vol. 62, no. 7, 2013, pp. 550–559.
62. I.G.Vladimirov, and I.R.Petersen, Coherent quantum filtering for physically realizable linear quantum plants, Proc. European Control Conference, IEEE, Zurich, Switzerland, 17-19 July 2013, pp. 2717–2723.
63. D.F.Walls, and G.J.Milburn, Quantum Optics, 2nd Ed., Springer, Berlin, 2008.
64. P.Whittle, Risk-sensitive linear/quadratic/Gaussian control, Adv. Appl. Probab., vol. 13, 1981, pp. 764–77.
65. R.M.Wilcox, Exponential operators and parameter differentiation in quantum physics, J. Math. Phys., vol. 8, no. 4, 1967, pp. 962–982.
66. H.M.Wiseman, and G.J.Milburn, Quantum measurement and control, Cambridge University Press, Cambridge.
67. N.Yamamoto, and L.Bouten, Quantum risk-sensitive estimation and robustness, IEEE Trans. Automat. Contr., vol. 54, no. 1, 2009, pp. 92–107.
68. K.Yosida, Functional Analysis, 6th Ed., Springer, Berlin, 1980.
69. G.Zhang, and M.R.James, Quantum feedback networks and control: a brief survey, Chinese Sci. Bull., vol. 57, no. 18, 2012, pp. 2200–2214,
arXiv:1201.6020v2 [quant-ph], 26 February 2012.

A Commutator of quadratic polynomials of quantum variables satisfying CCRs
Since the commutator adξ , associated with a given operator ξ , is a derivation, then repeated application of this property leads to
[ab, cd] = [ab, c]d + c[ab, d]
= ([a, c]b + a[b, c])d + c([a, d]b + a[b, d])
= [a, c]bd + a[b, c]d + c[a, d]b + ca[b, d]

(A1)

for any operators a, b, c, d. Hence, if these operators satisfy CCRs (that is, their pairwise commutators are identity operators up to scalar factors),
then the right-hand side of (A1) is a quadratic polynomial of a, b, c, d. Therefore, quadratic polynomials of such operators form a Lie algebra with
respect to the commutator. For the purposes of Section 4, we will provide a version of (A1) for quadratic polynomials of operators satisfying the
CCRs.
Lemma 4 Suppose a := (a j ), b := (bk ), c := (c` ), d := (dr ) are vectors of self-adjoint operators which satisfy the CCRs
 
 



a  T T
[a, cT ] [a, d T ]
Θ11 Θ12
, c d
:=
=
2i
,
b
Θ21 Θ22
[b, cT ] [b, d T ]

(A2)

where Θ jk are real matrices. Also, let F := ( f jk ) and G := (g`r ) be appropriately dimensioned complex matrices which specify the bilinear forms
aT Fb = ∑ f jk a j bk ,
j,k

cT Gd = ∑ g`r c` dr .

(A3)

`,r

Then their commutator is also a quadratic polynomial of the operators which is computed as

 
0
0
FΘ21 G
a
 T T T
T
T
T

0
0
F Θ11 G b .
[a Fb, c Gd] = 2i a b c
T F T GΘ T F
d
GΘ22
0
12
In the case a = b (when the vectors a and b are identical), (A4) takes the form

 


0
S(F)Θ11 G a
[aT Fa, cT Gd] = 4i aT cT
,
T
d
GΘ12 S(F)
0
where S(F) := 21 (F + F T ) denotes the symmetrizer of square matrices.

(A4)

(A5)


Proof By recalling the bilinearity of the commutator, applying (A1) to the operators a j , bk , c` , dr in (A3) and using the CCRs (A2), it follows that
[aT Fb, cT Gd] =

∑

f jk g`r [a j bk , c` dr ]

j,k,`,r

=

∑


f jk g`r [a j , c` ]bk dr + a j [bk , c` ]dr + c` [a j , dr ]bk + c` a j [bk , dr ]

j,k,`,r


T
T T
= 2i bT F TΘ11 Gd + aT FΘ21 Gd + cT GΘ12
Fb + cT GΘ22
F a .

(A6)

Quadratic-exponential Functionals and Large Deviations Estimates for Linear Quantum Systems

33

The right-hand side of (A6) is a quadratic function of the operators whose vector-matrix form is given by (A4). If a = b, then F in (A3) is a square
matrix, and the matrices in (A2) satisfy Θ1k = Θ2k for every k = 1, 2. In this case, (A6) reduces to

T
T T
[aT Fa, cT Gd] = 2i aT (F TΘ11 + FΘ21 )Gd + cT G(Θ12
F +Θ22
F )a

T
= 2i aT (F T + F)Θ11 Gd + cT GΘ12
(F + F T )a

T
= 4i aT S(F)Θ11 Gd + cT GΘ12
S(F)a ,
which establishes (A5).

Similar commutation relations hold for bilinear forms of annihilation and creation operators (see, for example, [25, Appendix B] and [54,
Lemma 4.2]) and are used in the context of Schwinger’s theorems on exponentials of such forms [9].

B Covariance of quadratic functions of Gaussian quantum variables
For the purposes of Section 5, we will need the following lemma on the covariance of bilinear forms of Gaussian quantum variables.
Lemma 5 Suppose the self-adjoint quantum variables, constituting the vectors a, b, c, d in Lemma 4, are in a zero-mean Gaussian state. Then the
covariance of the bilinear forms in (A3), specified by complex matrices F and G, can be computed as
E D
E
D
cov(aT Fb, cT Gd) = Fcov(b, d), cov(a, c)G + Fcov(b, c), cov(a, d)GT .

(B1)

In the case when a = b and c = d, the relation (B1) reduces to
D
E
cov(aT Fa, cT Gc) = 2 S(F), cov(a, c)S(G)cov(a, c)T ,

(B2)


where S is the symmetrizer.
Proof The mean value of the first bilinear form in (A3) is computed as
E(aT Fb) = ∑ f jk E(a j bk ) = Tr(F T E(abT )) = F, cov(a, b) ,

(B3)

j,k

where E(abT ) = cov(a, b) since the underlying quantum variables are assumed to have zero mean values. By a similar reasoning, the second bilinear
form in (A3) has the mean value
E(cT Gd) = G, cov(c, d) .

(B4)

Application of the Wick-Isserlis theorem [31, Theorem 1.28 on pp. 11–12] (see also [24] and [39, p. 122]) to the fourth-order mixed moment of the
zero-mean Gaussian quantum variables a j , bk , c` , dr leads to
E(aT FbcT Gd) =

∑

f jk g`r E(a j bk c` dr )

j,k,`,r

=

∑


f jk g`r cov(a j , bk )cov(c` , dr ) + cov(a j , c` )cov(bk , dr ) + cov(a j , dr )cov(bk , c` )

j,k,`,r

= Tr(F T cov(a, b))Tr(GT cov(c, d)) + Tr((Fcov(b, d))T cov(a, c)G) + Tr((Fcov(b, c))T cov(a, d)GT )
D
E D
E
= F, cov(a, b) G, cov(c, d) + Fcov(b, d), cov(a, c)G + Fcov(b, c), cov(a, d)GT .

(B5)

The covariance of the bilinear forms (A3) can now be computed by combining (B3)–(B5) as
cov(aT Fb, cT Gd) = E(aT FbcT Gd) − E(aT Fb)E(cT Gd)
D
E D
E
= Fcov(b, d), cov(a, c)G + Fcov(b, c), cov(a, d)GT ,
which establishes (B1). Application of (B1) to the particular case a = c and c = d leads to
D
E D
E
cov(aT Fa, cT Gc) = Fcov(a, c), cov(a, c)G + Fcov(a, c), cov(a, c)GT
D
E
= 2 Fcov(a, c), cov(a, c)S(G)
= 2 F, cov(a, c)S(G)cov(a, c)T
D
E
= 2 S(F), cov(a, c)S(G)cov(a, c)T ,

(B6)

thus proving (B2). In (B6) use has also been made of the symmetry of the matrix cov(a, c)S(G)cov(a, c)T , the orthogonality of the subspaces of
symmetric and antisymmetric matrices, and the fact that the symmetrizer commutes with the complex conjugation: S(F) = S(F).


34

Igor G. Vladimirov, Ian R. Petersen, Matthew R. James

C An averaging lemma with multifactor convolutions
For the purposes of the cumulant growth rate results of Theorem 4 in Section 6, we provide an averaging lemma which involves multifactor
convolutions (145).
Lemma 6 Suppose f1 , . . . , fr are appropriately dimensioned complex matrix-valued functions on the real line which are bounded and absolutely
integrable, with r > 2. Then
lim

t→+∞

r

1 Z

−
→

t

k=1

∏ fk (tk − tk+1 )dt1 × . . . × dtr
[0,t]r



= ( f1 ∗ . . . ∗ fr )(0)
1
2π

=
where tr+1 := t1 , and Fk (λ ) :=

R +∞
−∞

r
Z +∞ −
→

∏ Fk (λ )dλ ,

(C1)

−∞ k=1

fk (t)e−iλt dt denotes the Fourier transform of fk .



Proof Since the integrand on the left-hand side of (C1) depends on the integration variables t1 , . . . ,tr only through their differences, they can be
translated so as to represent the integral in the form
Z

r

−
→

Z t

k=1

0

∏ fk (tk − tk+1 )dt1 × . . . × dtr =
[0,t]r

gt (t1 )dt1 ,

(C2)

where
r

−
→

Z

∏ fk (tk − tk+1 )dt2 × . . . × dtr
[0,t]r−1

gt (t1 ) :=

k=1

r

−
→

Z

=

∏ fk (τk − τk+1 )dτ2 × . . . × dτr .
[−t ,t−t ]r−1
1

1

(C3)

k=1

Here, use is made of the new integration variables τk := tk − t1 for all k = 2, . . . , r together with the corresponding convention τ1 = τr+1 = 0.
The right-hand side of (C3) is organized as the convolution ( f1 ∗ . . . ∗ fr )(0) (evaluated at 0) except that the integration is restricted to the cube
[−t1 ,t −t1 ]r−1 . Now, the fulfillment of the inclusions t1 ∈ [0,t] and (τ2 , . . . , τr ) ∈ [−t1 ,t −t1 ]r−1 is equivalent to t1 belonging to the (possibly empty)
interval
− min(0, τ2 , . . . , τr ) 6 t1 6 t − max(0, τ2 , . . . , τr )
(C4)
of length

max 0, t − max(0, τ2 , . . . , τr ) + min(0, τ2 , . . . , τr ) = tht (τ2 , . . . , τr ).

(C5)


ht (τ2 , . . . , τr ) := χt max(0, τ2 , . . . , τr ) − min(0, τ2 , . . . , τr )

(C6)

Here,
inherits from the function χt in (105) the properties of being bounded by and convergent to 1 as t → +∞ for any given τ2 , . . . , τr ∈ R. By combining
(C2)–(C6), it follows that
1
t

Z

r

−
→

∏ fk (tk − tk+1 )dt1 × . . . × dtr =
[0,t]r

r

Z

k=1

→

Rr−1

Z

−
→
ht (τ2 , . . . , τr ) ∏ fk (τk − τk+1 )dτ2 × . . . × dτr
k=1

r

−
→

∏ fk (τk − τk+1 )dτ2 × . . . × dτr ,

Rr−1 k=1

as t → +∞,

(C7)

where Lebesgue’s dominated convergence theorem is applicable since the functions f1 , . . . , fr are bounded and absolutely integrable (and hence, so
are their convolutions). The limit in (C7) is ( f1 ∗ . . . ∗ fr )(0), which establishes the first of the equalities (C1), with the second of them following
from the convolution theorem.


