Highly accurate Gaussian process tomography with geometrical sets of coherent states

We propose a practical strategy for choosing sets of input coherent states that are near-optimal for reconstructing single-mode Gaussian quantum processes with output-state heterodyne measurements. We first derive analytical expressions for the mean squared-error that quantifies the reconstruction accuracy for general process tomography and large data. Using such expressions, upon relaxing the trace-preserving constraint, we introduce an error-reducing set of input coherent states that is independent of the measurement data or the unknown true process -- the geometrical set. We numerically show that process reconstruction from such input coherent states is nearly as accurate as that from the best possible set of coherent states chosen with the complete knowledge about the process. This allows us to efficiently characterize Gaussian processes even with reasonably low-energy coherent states. We numerically observe that the geometrical strategy without trace preservation beats all nonadaptive strategies for arbitrary trace-preserving Gaussian processes of typical parameter ranges so long as the displacement components are not too large.

the best TP strategy so long as the process displacement components are not very large.
After some background introduction to the general formalism of Gaussian processes in Sec. II, Sec. III shall be devoted to the explanation and derivation of the MSE formulas for both TP and non-TP reconstruction methods. With the aid of these formulas, Sec. IV then proceeds with the construction of geometrical input states. Finally, Sec. V compares the geometrical strategy with existing common nonadaptive input-state strategies for realistic CPTP Gaussian processes.

II. CHARACTERIZATION OF GAUSSIAN PROCESSES
A physical quantum process Φ transforms an input state ρ IN into the output state ρ OUT = Φ[ρ IN ].
A standard operational description for the quantum process Φ makes use of the Choi-Jamiołkowski formalism, which essentially states that all information about Φ is encoded into a positive operator (ρ Φ ). Additionally, we say that Φ is Gaussian if it possesses a two-mode Gaussian quasidistribution. In this case, it is convenient to represent Φ by its Husimi Q function and a real constant c 0 . Here Z Z Z = (z z z 1 z z z 2 ) consolidates the complex variables labeling the process input [z z z 1 = (z 1 z * 1 ) ] and output [z z z 2 = (z 2 z * 2 ) ] modes. Gaussian process tomography pertains to the characterization of any given unknown ρ Φ on the premise that its Q Φ is Gaussian. The connection between ρ Φ and Q Φ is made by heterodyne measurements [29][30][31][32][33][34][35][36] that sample the overcomplete set of coherent states {|z z|} to probe the output state ρ OUT = tr 1 {ρ IN ⊗ 1 ρ Φ }, where the transposition is defined for the Fock basis in which all matrices are written in this article. Apart from directly recovering the Q-function parameters, these measurements are also known to give a smaller MSE for characterizing covariance matrices of Gaussian and broad classes of non-Gaussian quantum states compared to its homodyne counterpart [37][38][39]. Another reason for this choice of measurements is that when coherent input states {|α α|} are used, the heterodyne measurement is equivalent to a direct sampling of the process Q function, as Q Φ (α, α * , z, z * ) = α|ρ OUT |α .
For a completely-positive (CP) Φ (ρ Φ ≥ 0), A A A, B B B and c 0 are constrained such that Q Φ is positive and square-integrable. One way to identify these constraints systematically is by reverting to the real phase-space representation: . This is done by recognizing that the transformations Z Z Z = U U U R R R, A A A = U U U † A A AU U U and B B B = U U U † B B B are exacted with the unitary matrix U U U = 1 1 1 ⊗U U U 0 and U U U 0 =   1 i

It is now clear that the conditions
A A A ≥ 0 and A A A ≥ 0 are necessary for Q Φ to be real and square-integrable. These give a total of 15 independent real parameters, that is 10 from A A A, 4 from B B B, and c 0 . We may parametrize A A A and B B B as A useful Φ in quantum information theory is typically also TP (tr{ρ IN } = 1 = tr{ρ OUT }). Under this constraint, for an invertible A A A 3 , it is shown in Appendix A that 6 of the 15 real parameters are fixed by the rest inasmuch as As a simple example, if we consider a beam splitter that transforms a pair of input mode operators a and b into the pair of output operators c = a cos θ + b sin θ and d = −a sin θ + b cos θ , then the relevant Choi-Jamiołkowski operator for a single output mode (c) clearly describes a CPTP process and possesses the Q function Q Φ = exp(−|z 1 | 2 (cos θ ) 2 − |z 2 | 2 + z 1 z 2 cos θ + z * 1 z * 2 cos θ ) [45]. In this case, a consistency check gives 0 0 0 and c 0 = 0, with σ σ σ x being the usual Pauli x matrix in the standard basis.
If the TP constraint is absent from the process reconstruction, the parameter c 0 is not estimable since any experimental data can only recover ρ Φ uniquely up to a constant multiple [46,47], such that Φ may only be fully characterized up to its operator trace. Therefore, the complete characterization of a general single-mode Gaussian Φ requires 14 real recoverable parameters: x x x = (a 1 a 2 b 1,r b 1,i b 2,r b 2,i c 1,r c 1,i c 2,r c 2,i g 1,r g 1,i g 2,r g 2,i ) , (2.4) where the subscripts r and i denote the real and imaginary parts of a complex parameter. Formally, the coherent-state sampling measurements by heterodyning gather raw data sampled from the Q function Q OUT (x x x; z 2 ≡ z, z * 2 ≡ z * ) of ρ OUT (originating from a given ρ IN ) that encodes x x x. Numerical techniques are then used to obtain the estimator x x x (distinguished from the true parameter by a caret).
In this work, we focus on studying the accuracy of x x x for a given unknown x x x. This may be quantified by the MSE ( x x x −x x x) 2 , where the overline denotes an average over all possible data of a fixed total sample size. For an analytical study, we shall investigate the asymptotic expression of the MSE that applies to typical tomography situations involving large datasets.

A. Relaxation of the TP constraint
We first investigate the case where coherent states are used as input states for characterizing an unknown, generally non-TP Gaussian Φ. Given a ρ IN = |α α|, the output Gaussian Qfunction Q OUT can be written in the form where x x x = (x x x c 0 ) carries the non-estimable c 0 . In Gaussian process tomography, x x x may be extracted from a discretized system of equations governed by (3.1). The latter is established by first sending J input coherent states {|α j α j |} J j=1 , next performing heterodyne measurements on all corresponding output states ρ OUT (x x x; z k , z * k ) up to proper normalization as a consequence of binning.
For LI to be successful, the system u u u = V V Vx x x must be informationally complete (IC), that is, This implies that measurement data collected with such a set of input states uniquely characterize the unknown Gaussian process.
In equivalent linear-algebraic terms, an IC set of linearly independent input states gives rise to To understand why this is the case, we observe that as 6 out of the 15 terms in v v v do not depend on z, when J < 6, there naturally exists at least one null right eigenvector e e e for V V V of the form e e e = (e 1 0 e 2 e 3 0 0 e 4 e 5 0 0 0 0 0 0 e 6 ) , where the 6-dimensional (e 1 e 2 e 3 e 4 e 5 e 6 ) is orthogonal to (|α j | 2 2α j,r 2α j,i 2(α 2 j ) r 2(α 2 j ) i 1) for any amplitude α j . This observation is therefore consistent with the alternative arguments in [45].
In realistic scenarios, the log-probability column u u u is to be replaced by the column of relative log-frequencies u u u = (− log ν jk ) that reflects the physical relative photodetection counts. Note that ∑ k ν jk = 1, and ν jk → p jk as N 1 in a statistically consistent setting. As a consequence, the LI procedure that now handles these noisy data ν jk should be modified. As the counts are noisy with statistical fluctuation, for any finite number of sampling copies N per input state, there very likely exist entries in u u u that are infinite (ν jk = 0 for some j and k), especially when the corresponding Q-function magnitudes are small. To cope with statistical noise in LI, one may consider only finite entries of u u u. After some statistical reasoning (see Appendix B), we obtain the asymptotic where we note that ∑ k p jk = 1 are the normalized true probabilities related to p jk through p jk = p jk / ∑ k p jk . Figure 2 illustrates the positive match between the theoretical expression in (3.2) and simulation results for a given Gaussian process. The real and imaginary parts of the complex input coherent-state amplitude α j = α j,r + i α j,i are chosen from the closed interval [−L, L], where the influence of L on the characterization quality of Φ is explored. We assume that γ j ≡ ∑ k p jk are known with small statistical fluctuation up to a scalar multiple.
We stress that the LI estimator x introduced here, while useful as a formalism for an analytical grasp of the actual characterization problem, usually does not lead to a physical process estimator Φ, since the inversion procedure pays no attention to the positivity requirement for the estimated complex A A A matrix. Numerically, it is possible to enforce such a positivity constraint in LI, in which case the resulting estimator will have some statistical bias and an MSE that deviates slightly from the expression in (3.2). One may also choose to perform LI followed by a projection onto the real and positive A A A -space in the real phase-space representation, as previously discussed in Sec. II, to obtain a sufficiently good physical estimator x. Supposing that the estimated real matrix While LI with positivity constraint and the projection method give highly similar estimators for sufficiently large N, the scaling behaviors in JN for both methods generally vary. To put things on firmer statistical grounds, more meaningful estimators, such as the maximum-likelihood (ML) estimators, should be considered. In the context of non-TP process characterization, ML asymptotically gives very similar reconstructions to LI under the physical process constraints. Section III B provides an explanation regarding this connection and presents a recipe for the ML reconstruction prescription.

B. Imposition of the TP constraint
If the unknown Gaussian process Φ is TP, then the constraints specified in (2.3) dictate that 9 parameters are enough to characterize Φ: We note that in this case, A A A ≥ 0 so long as A A A 3 > 0, since we may write with well-defined matrix square-roots. This also implies that A A A is rank-2 and that a 2 2 − 4|c 2 | 2 > 0 is the only necessary and sufficient positivity condition for a CPTP Gaussian Φ as no other constraints are imposed on A A A 2 and b b b 2 .
Because of the nonlinear dependence on the 9 parameters in the exponent of the output Q function in accordance with (2.3), LI is no longer applicable as it only works with a highly specific form of the output Q function stated in (3.1). Instead, statistical method is usually a more favorable option to infer x x x from the collected data. A popular method is to maximize the log-likelihood function log L = ∑ jk ν jk log(p jk / ∑ j k p j k ) that takes the multinomial form when each output state is measured with N sampling copies of heterodyne detection independently, subject to the positivity constraint of A A A 3 ≥ 0 0 0. The log-likelihood log L may in general be a nonconvex function of the TP Gaussian-process parameters, so standard numerical techniques might be needed to search for its global maximum for optimal accuracy.
We emphasize that the ML scheme may be applied to any tomographic situation, which evidently includes the characterization of non-TP processes. In this context, with respect to the where we again assume that the γ j s can be determined through calibration procedures up to a multiplicative constant and are not part of the statistical consideration. Maximizing log L involves for the parameter x x x and µ = ∑ j γ j . If the solution to µ p jk / ∑ j k p j k = ν jk exists under the physical constraints of the parameter estimator x x x for which the corresponding process estimator The asymptotic MSE expression for the maximum-likelihood (ML) estimator x x x ML for x x x may be approximately calculated by assuming that N 1 per output state is large enough so that JK × 9 matrix V V V TP such that its ( j, k)th row is equal to the 1 × 9 row where the definitions of all auxiliary matrices are given in Appendix B. The principle of small variations thus states that V V V TP δx x x = δu u u. It follows that the asymptotic MSE for estimating the 9 independent CPTP parameters using the constrained ML method is approximately where Y Y Y is as specified in Eq. (3.2).
The ML estimator x x x ML is to be constrained by the positivity of A A A 3 and is typically a biased estimator. In general, there could still be a small difference between the ML MSE and the righthand side of Eq. (3.8) that is obtained from operator derivatives that assume the existence of open sets without parameter boundary constraints. Barring this technical issue, Fig. 3 shows that (3.8) can still serve as a pretty good estimate for the actual MSE.

IV. GEOMETRICAL SET OF INPUT COHERENT STATES
The LI procedure discussed in Sec. III A hinges on the existence of V V V − . This is again synonymous to having a V V V with no null right eigenvectors, thereby ruling out sets of coherent states with identical amplitudes |α j | ≡ |α| as candidates for LI, since they result in at least one such null eigenvector, namely e e e ∝ (−1/|α| 2 0 . . . 0 1) , for any J. Therefore, choices that include the set of coherent states with complex amplitudes that form a ring of radius r in phase spaceα j = r e 2πi j/J for r > 0 and 1 ≤ j ≤ J-are non-IC and shall result in failures of the LI scheme.
In the regime of ML estimation, these symmetric sets of coherent states give a convex set of estimated parameters that are consistent with the ML probabilities obtained from the measurement data. One therefore cannot obtain a unique parameter reconstruction with such input states.
A general observation from Figs. 2 and 3 is that input coherent states with lower phase-space energies (L = 1 for instance) also tend to give larger average MSE values. This is because a set of low-energy coherent states are typically quite closely packed in phase space and their Gram matrix G G G can at times be ill-conditioned, that is, its smallest eigenvalue can be very close to zero. Colloquially, low-energy coherent states are not very linearly independent. However, with appropriate optimization strategies, low-energy input coherent states can still be highly effective in characterizing any unknown Gaussian process with significantly higher accuracy than using random low-energy input states, thereby allowing us to maximally utilize these energy-efficient resources.
Regardless of whether or not the TP constraint is imposed when reconstructing Φ, the best performance of any Gaussian-process characterization is ultimately tied to the optimal value of the figure of merit in question. For our case, this is quantified by the minimum of the MSE. The truly optimal set of input states that minimizes the MSE according to either Eq. (3.2) or (3.8) requires the knowledge of the unknown Gaussian process. Such a set is therefore operationally unobtainable.
We introduce a solution to approximately minimize the MSE without such knowledge by first noting that since Y Y Y ≥ 0, a variant of the Cauchy-Schwarz inequality for positive matrices reads For low-energy applications, L = 1 is sufficient for precise Gaussian-process tomography.
where the approximation in the second line of the calculation is valid for sufficiently large M, so that p jk 1 and Y jk, j k ≈ N 2 p jk δ j, j δ k,k . Under this approximation, it is clear that the upper bound of the asymptotic MSE is independent of Φ and is therefore a purely geometrical term that solely depends on the collective phase-space arrangement of the input coherent states. Therefore, such geometrical sets of input states that minimize the MSE Cauchy-Schwarz upper bound can be universally defined for any Gaussian process since they are independent of the measurement data and the unknown process. Figure 4 illustrates some desirable properties of geometrical input states. In particular, when J = 6 (minimal case), the geometrical set is unique up to a collective rotation, whereas arrangements can vary for J > 6, at times with the possibility of two or more coherent states being very close to each other owing to overcomplete redundancy. The recipe for deriving such geometrical sets of input states does not apply to the TP version of the asymptotic MSE in Eq. (3.8), since the TP constraint tangles the 9 independent Gaussian parameters in a highly nonlinear way that cannot be cleanly separated from the phase-space variables (V V V TP depends on these parameters).
As we shall see in Sec. V, for heterodyne detection, such a geometrical set of coherent states on average gives a nearly-optimal MSE when the TP constraint is lifted. Moreover, these special input states can also beat optimal input states that minimizes the MSE when trace preservation is imposed for certain classes of CPTP Gaussian processes.

For a given CPTP Gaussian process, both its Q-function parameters A A A and B B B are related to another set of parameters (refer to Appendix C) according to the maps
where U U U is the unitary matrix defined in Sec. II. The matrices X X X and Y Y Y effect the general transformations µ µ µ → X X Xµ µ µ + µ µ µ 0 and Σ Σ Σ W → Σ Σ Σ W = X X X Σ Σ Σ W X X X +Y Y Y on the mean (µ µ µ) and covariance (Σ Σ Σ W ) of the Wigner function describing an input Gaussian state, and The matrices X X X and Y Y Y may then alternatively be understood as functions of a complete set of relevant operations, namely phase shift (φ ), displacement (x 0 , p 0 ), squeezing (r, θ ), losses (χ < 1) and amplifications (χ > 1), and couplings to a Gaussian reservoir (n T , a T , θ T ): transmission rates above 0.9 [50]. Amplification can be a result of nonlinear processes, but more often, it arises as a consequence of feed-forward with non-unit gain [51,52]. Like displacement, it can in principle be arbitrary, but since strong amplification necessitates noise levels that destroys non-classical features of quantum states [53], it is usually kept low. As a conservative choice, we pick χ ∈ (0.1, 1.5) to run all numerical experiments. The last three terms collectively characterize the added noise, which has purely detrimental effect-adding one unit of vacuum noise to both quadratures is generally sufficient to extinguish any quantum properties of the state. The asymmetry coefficient a T ∈ [−1, 1] and phase θ T ∈ [0, π/2] of the added noise span the practical ranges.
The added-noise coupling n T ∈ [0, 1] was chosen such that it describes both quantum and classical channels.
Three groups of Gaussian CPTP processes are considered here, whose physical parameters are tabulated in Tab 5(a,b,c) and Figs. 6(a,b,c)], GML is the optimal choice for efficient process-parameter reconstruction with low-energy coherent states (L = 1).
Interesting dynamics reveal themselves for the completely random CPTP processes in Gp. 3, where we find that, again, GML beats all strategies when all parameters are arbitrarily chosen with the displacement ranges −2 ≤ x 0 , p 0 ≤ 2 obeyed. For the tested random processes, Figs. 5(c) and 6(c) highlight rather comparable performances between GML and the optimal BML (TP). How-ever, as shown in Fig. 7, if we enlarge the displacement ranges, we find that the best TP input-state strategy can outperform the rest. On the other hand, when these ranges are reduced, GML reconstructs process parameters with much better accuracies than for default ranges as in Figs. 5(c) and 6(c). This leads us to conjecture that the additional reconstruction bias introduced by the TP constraint, on top of that from the CP constraint, apparently has beneficial merits for more correctly singling out estimators that are near the true CPTP Gaussian processes that perform large displacing operations; whereas processes with stronger second-moment manipulating features relative to first-moment displacements are still better characterized with non-TP input-state strategies as they appear to be more robust against noise, in which case GML is the optimal choice.

VI. CONCLUSION
There have been many studies related to achieving the quantum limits of parameter estimation.
We have taken a different route instead and investigated an experimentally feasible and tomographically efficient way to characterize Gaussian quantum processes. Using heterodyne measurements and input coherent states, we introduced a simple strategy of constructing geometrical sets of input coherent states that effectively optimizes the process parameters' mean squared-error. These geometrical input states are demonstrated to outperform the best nonadaptive input-state strategy in terms of the mean squared-error for typical CPTP processes that do not carry out large displacement operations. This permits us to utilize coherent states of low energies as sufficient convenient resources to achieve very low mean squared-errors in the reconstructed process parameters.
We also observe that if the unknown Gaussian process has large displacing features on input states, input-state strategies that imposes the trace-preserving constraint, on average, give more accurate process estimators than the geometrical strategy where this constraint is relaxed.
In the course of acquiring these results, we also obtained asymptotic analytical expressions for the process-parameter mean squared error that were previously not discussed in the quantum process The TP constraint imposed on a quantum process Φ is defined by the partial-trace relation tr 2 {ρ Φ } = 1 for its Choi-Jamiołkowski operator ρ Φ . Using the overcompleteness property of the coherent states, (dβ )|β β |/π = 1, this partial-trace relation is expressed as where a a a = (a a † ) is the column of ladder operators. The normal-order expression in (A1) may be unraveled with two other overcomplete sets of coherent states, The Gaussian integration in β can be evaluated by using the general formula for Re{a} ≥ |c 1 + c * 2 | (see [44,54]). At times, it is much easier to cope with the matrix form of the integration: After an application of (A4) and some reorganization of terms, we have |α α|α e − α α αW W W α α α +y y y † α α α + w 0 α | , It therefore follows that the sufficient conditions for Φ to be TP are W W W = 0 0 0, y y y = 0 0 0 and w 0 = 0.
To show that these are necessary conditions, we rewrite (A5) as : e −a a a † W W Wa a a +y y y † a a a + w 0 : = 1 (A6) and make use of Eqs. (1.11)-(1.13) in [55] to convert the normal-order form on the left-hand side of (A6) into : e −a a a † W W Wa a a +y y y † a a a : = 1 κ e −a a a † W W WR R R 1 T T T −1 R R R 1 a a a +y y y † R R R 1 T T T −1 R R R 1 a a a , where κ = det sinc(σ σ σ y ξ ξ ξ ) −1/2 det{T T T } 1/2 exp(η η η (T −1 1 1 T −1 1 1 T −1 1 1)ξ ξ ξ −1 η η η) , It is now obvious that the earlier sufficient conditions are also necessary for the Gaussian form on the right-hand side of (A7) to be the identity operator, as R R R 1 T T T −1 1 R R R 1 is non-singular.
as p jk → 0. The final expressions stated in (3.2) hence follow suit.
The derivation of the MSE expression for the ML reconstruction scheme is only slightly more technical. In calculating the increment δ log L , variations in the respective parameter matrices and column lead to where the gradient components M M M 1, jk , M M M 2, jk and M M M 3, jk can be derived using the simple variational At this stage, we introduce the vectorization notation |Y Y Y that refers to the column formed by stacking all columns of Y Y Y in the computational-basis representation, and its dual |Y Y Y † ≡ Y Y Y † |. These higherdimensional objects relate to the trace inner product of two matrices Next, we identify the essential transformations to calculate the 9 × 1 δx x x: All these auxiliary matrices form the components needed to establish Eq. (3.8).

Appendix C: Physical features of CPTP Gaussian processes
In phase-space representation, the second-moment matrix A A A of the process Q function exp (−R R R A A A R R R+ B B B R R R + c 0 ) is related to the Q-function covariance Σ Σ Σ Q as A A A = Σ Σ Σ −1 Q /2 by definition. Since the Q function is a Gaussian convolution of the Wigner function, this Q-function covariance can in turn be obtained from the Wigner-function covariance Σ Σ Σ W in accordance with Σ Σ Σ Q = Σ Σ Σ W +1 1 1/2. We may equivalently regard Σ Σ Σ W as the output of a general covariance transformation on the otherwise idle channel covariance Σ Σ Σ 0 : where X X X and Y Y Y precisely effect such a transformation on the covariance matrix Σ Σ Σ defined for the Gaussian Wigner function of the input state inasmuch as Σ Σ Σ W = X X X Σ Σ Σ W X X X +Y Y Y . The matrices X X X and Y Y Y are functions of the operational parameters defined in Sec. V. Although their functional forms may not be uniquely established, the transformed covariance Σ Σ Σ W must satisfy the inequality Σ Σ Σ W + iΩ Ω Ω/2 ≥ 0 0 0 for Ω Ω Ω = which preserves Heisenberg's uncertainty relations [11]. To confirm that the decomposition in (5.3) for X X X and Y Y Y precisely obeys this inequality, recall that X X X = χ S S SR R R is proportionally symplectic-X X X Ω Ω ΩX X X = χ 2 R R R S S S Ω Ω ΩS S SR R R = χ 2 Ω Ω Ω. On the other hand, Y Y Y = |1 − χ 2 |1 1 1/2 + Y Y Y 0 is a linear combination of a symplectic multiple and a positive matrix Y Y Y 0 . Thus, Σ Σ Σ W + iΩ Ω Ω/2 = |1 − χ 2 |1 1 1/2 − i(χ 2 − 1)Ω Ω Ω/2 +X X X (Σ Σ Σ + iΩ Ω Ω/2)X X X +Y Y Y 0 , where the final two matrices are clearly positive. The first two matrices combine to give a projector multiple for any χ. In terms of the standard Pauli matrix σ σ σ y , they are (1 − χ 2 )(1 − σ σ σ y )/2 and (χ 2 − 1)(1 +σ σ σ y )/2 for the respective ranges χ < 1 and χ > 1.
which is rank-2 due to the CPTP character (recall Sec. III B). Moreover, A A A 0 is a rank-2 projector- There is therefore strictly no valid Σ Σ Σ 0 for the idle process as A A A 0 is not invertible in the usual sense. The validity of (C1) would demand an understanding that A A A is truly invertible as soon as the TP constraint is violated, and that all parameter variations in A A A should encounter no discontinuities when transiting between non-TP and TP subspaces. Upon recognizing that the Choi-Jamiołkowski operator for the idle channel is represented by the superposition ∑ ∞ n=0 |nn , we may consider a more physical realization using a two-mode squeezed vacuum state of real finite squeezing strength ω 1: |nn (tanh ω) n+n n n | .
It follows that the relevant ω-deformed Q function spells Q Φ (ω) = α, z|Φ idle ω |α, z = e −|α| 2 − |z| 2 + (αz + α * z * ) tanh ω , from whence one reads off the full-rank matrices Obviously, in the limit ω → ∞, we revert to A A A → A A A 0 . Yet, any finite ω leads to a well-defined Σ Σ Σ t = (2A A A ) −1 | t=2ω − 1 1 1/2 in Eq. (5.2) after the substitution t = 2ω. That the inverse of A A A 0 for the idle process is ill-defined should now be even clearer from lim t→∞ Σ Σ Σ t = Σ Σ Σ 0 . Despite this divergence, the limit expression in (5.1) must tend to a finite matrix A A A, as to be expected from the physical point of view.