Inverse counting statistics for stochastic and open quantum systems: the characteristic polynomial approach

We consider stochastic and open quantum systems with a finite number of states, where a stochastic transition between two specific states is monitored by a detector. The long-time counting statistics of the observed realizations of the transition, parametrized by cumulants, is the only available information about the system. We present an analytical method for reconstructing generators of the time evolution of the system compatible with the observations. The practicality of the reconstruction method is demonstrated by the examples of a laser-driven atom and the kinetics of enzyme-catalyzed reactions. Moreover, we propose cumulant-based criteria for testing the non-classicality and non-Markovianity of the time evolution, and lower bounds for the system dimension. Our analytical results rely on the close connection between the cumulants of the counting statistics and the characteristic polynomial of the generator, which takes the role of the cumulant generating function.


Introduction
When probing a physical system we often face the problem that only a small part of its evolution is accessible to direct observation. Notably, in some cases, the entire available information consists of observable time-discrete events, indicating that a stochastic transition between two states of the system has been realized. Examples of such events are photons emitted by fluorescent ions [1] or nitrogen-vacancy centers [2], single electrons passing through quantum dots [3][4][5], individual steps of processive motor proteins [6] or product molecules generated in enzyme-catalyzed reactions [7][8][9]. Common to all is that the observation of the system is restricted to a simple counting process of certain discrete events by a detector.
In this paper we ask the question: What can be learned about the physical system if the counting statistics of these events is the only available information? To be specific, the longtime counting statistics is the discrete probability distribution p(n), expressing the probability that exactly n events occur within a sufficiently long time interval [0, t]. The distribution p(n) is conveniently parametrized in terms of its cumulants, where non-zero higher order cumulants indicate deviations from Gaussian behaviour. Our goal is to exploit these deviations in order to recover the hidden structure and evolution of the system from the observational data. Figure 1. A detector monitors a stochastic transition between two specific states (red arrow) of a system with N states connected by stochastic and quantum transitions (arrows and dashed lines). The long-time counting statistics of the events p(n) obtained from the detector signal D(t) contains sufficient information to reconstruct the generators governing the time evolution of the system. In answering the above question, we present several results that elucidate the relation between the evolution of the system and the observed counting statistics. Most importantly, we identify the characteristic polynomial of the time-independent generator, governing the time evolution, as the central object of our study. The characteristic polynomial is shown to have two essential properties: First, it can be reconstructed from a finite number of cumulants, and second, it takes the role of the cumulant generating function. As an immediate consequence, we find that for an open quantum system with dimension N at most 2(N 2 − 1) cumulants of the observed counting statistics are independent. Building on these results, we further design cumulant-based tests to identify non-classical and non-Markovian time evolutions, and to estimate the system dimension N. Moreover, we suggest an analytical method for reconstructing generators compatible with the observations. We find that the result of the reconstruction process is not necessarily unique, i.e., different generators may lead to identical observed counting statistics. By analogy to problems such as inverse scattering theory [10] we refer to our methods collectively as inverse counting statistics (ICS).
The theory of ICS is a new member of a family of methods to discern the properties of a system from observational data. Part of these methods have the objective of reconstructing specific system properties, such as quantum state tomography [11][12][13][14], which addresses the problem of efficiently determining the state of a physical system by measuring a complete set of observables. Various methods, on the other hand, are based on inequalities for selected system properties; these inequalities are tested against experimental observations. The most established test for the non-classicality of the evolution of a system is the Leggett-Garg inequality [15,16], whose experimental violation was demonstrated with nitrogen-vacancy defects [17] and which may be used to evidence non-classical electron transport in nanostructures [18]. The advantages of cumulant-based inequalities for probing non-classicality of current correlations in mesoscopic junctions and a related cumulant-based Bell test have been recently discussed in Refs. [19,20]. Furthermore, identifying non-Markovian evolutions has become a subject of increasing interest and several measures of non-Markovianity in quantum systems have been proposed [21][22][23]. For determining the dimension of quantum systems, the concept of dimension witness has been introduced [24] and experimentally applied [25]. The dimensionality of stochastic systems in the context of enzyme dynamics was shown to be directly related to the so-called randomness parameter [6,8,9].
The theory of ICS is of course related to the methods of full counting statistics (FCS) [26][27][28][29][30][31][32] and large deviation theory [33][34][35][36]. These methods have been successfully applied to a wide range of stochastic and quantum systems, notably in mesoscopic physics [32]. In fact, measurements of higher order cumulants of electronic current fluctuations have been achieved in tunnel junctions [37,38] and in Coulomb blockaded quantum dots [3][4][5]. In general, FCS and large deviation theory are aimed at determining the counting statistics of observable events if the generator of the time evolution and all system parameters are known a priori. Therefore these methods are not directly applicable to the inverse problem considered here.
The structure of this paper is as follows: In section 2 we introduce the underlying model and recall the basics of the theory of FCS. In section 3 we give the details of the methods of ICS by establishing the connection between the characteristic polynomial and the cumulants of the counting statistics. Furthermore, we describe the cumulant-based tests and the analytical method for reconstructing the generator of the time evolution. To illustrate the capacities of ICS we present specific reconstruction examples in section 4. We end with the conclusions in section 5.

Model and prerequisites
While keeping the discussion broad we consider a specific system for clarity. This system consists of N orthogonal states, connected by stochastic and quantum transitions (cf. Fig 1). The state of the system is described by the density operator ρ(t) whose time evolution is governed by the master equation with L the time-independent Markovian generator. The formal solution of Eq. (1) is given by ρ(t) = e Lt ρ(0). We focus specifically on the long-time limit of the evolution of the system and presume that Eq. (1) has a unique steady state ρ s defined by Lρ s = 0, i.e., the generator L has a unique zero eigenvalue. We want our model to cover two scenarios: In the classical limit, the system evolves according to a time-continuous Markov process [39] with only stochastic transitions, whereas in the general case we deal with an open quantum system [40,41]. The model can also be seen as a continuous-time extension of hidden quantum Markov models studied in Ref. [42]. The generator L in Lindblad form acts on ρ as where L ij and L † ij are Lindblad operators and { , } stands here for the anticommutator. The Lindblad operators are of the form L ij = √ κ ij |i j| and describe stochastic transitions from state |j to state |i with time-independent rates κ ij ≥ 0. The time-independent Hamiltonian operator H is expressed as a Hermitian matrix in the basis {|i }, fixed by the stochastic transitions. A priori, this model has 2N 2 − N real-valued parameters: N 2 − N transition rates κ ij , (N 2 − N)/2 complex-valued off-diagonal elements of H, and N real-valued diagonal elements of H. For many systems of interest, part of the stochastic or quantum transitions between states are identically zero; the non-zero transitions are conveniently depicted as a graph (cf. Fig 1). We assume that the detector monitors a single stochastic transition |j * → |i * between two selected states |j * and |i * with perfect efficiency. The detector produces a timecontinuous signal D(t), containing the number n ≥ 0 of events observed up to the time t (cf. Fig. 1). The combined state of the system and detector after the n-th incoherent event is given by the n-resolved density operator ρ n (t), where the additional index specifies the state of the detector [43].
The presented model is applicable to many interesting settings, in particular to stochastic and quantum mechanical transport problems [32,39]. In this case, N − 1 states describe different particle configurations and the remaining state is identified with the empty state of the system. Transitions that involve a change in the number of particles are interpreted as particle transfers between the system and external reservoirs, and the detector measures particle currents and fluctuations.

Established results from FCS
Let us briefly recall useful results from the theory of FCS (see Ref. [31] for details) and define quantities of interest for the rest of the paper. The density operator of the system and detector ρ n (t) satisfies the n-resolved master equation ∂ t ρ n (t) = L (0) ρ n (t) + L (1) ρ n−1 (t), where the generator has been separated into two parts as L = L (0) + L (1) . In terms of matrix elements, L (1) = L i * j * ρL † i * j * contains the only off-diagonal element corresponding to the transition |j * → |i * and L (0) the remaining elements. The finite-difference equation for ρ n (t) is solved by using the discrete Laplace transform ρ ξ (t) = n ρ n (t)e ξn which obeys the equation where the deformed generator L ξ = L (0) + e ξ L (1) has been introduced. The Fourier transform is often used instead of the Laplace transform, which corresponds to the replacement ξ → iξ throughout the paper. A crucial observation of FCS and large deviation theory is that the cumulant generating function G(ξ) = log E e ξD , where E · stands for the expectation, is connected to the eigenvalue of smallest magnitude λ(ξ) of the deformed generator L ξ through the relation lim t→∞ G(ξ)/t = λ(ξ) [33,36]. In other words, λ(ξ) corresponds to the unique zero eigenvalue of L in the limit ξ → 0. In the long-time limit, the cumulants C ν of the probability distribution p(n) of the detector variable D(t) are then given by Using Eq. (4) therefore allows us to determine the counting statistics of the detector variable D(t), in particular its average E D = C 1 and variance Var(D) = C 2 , provided that the deformed generator L ξ is known. Since all cumulants C ν increase linearly with time t, we introduce the scaled cumulants c ν = C ν /t. Corresponding results for the purely classical Markov process are obtained by restricting the previous derivations to the occupation probabilities ρ ii . It is then convenient to introduce the probability vector p = {ρ ii , i = 1, . . . , N} and the stochastic generator L cl acting on p. In practice, L cl and the quantum generator L qm are represented by N × N and N 2 × N 2 matrices, respectively, and ρ = {ρ ij , i, j = 1, . . . , N} is a vector in Liouville space with N 2 components. Unless necessary, we will not explicitly distinguish between stochastic and open quantum systems and consider a generic generator of dimension M × M.

Inverse counting statistics
We now present the general theory of ICS, which enables us to analyze the system by using the information provided by measured cumulants. The main results of ICS are of two kinds: The first ( §3.3) is formulated as a cumulant-based test for non-classicality, non-Markovianity and system dimension N; the second ( §3.4) aims at the reconstruction of generator of the time evolution. Both are based on the close relation between the characteristic polynomial of the generator and the cumulants, which will be established first.

Characteristic polynomial of the generator
Most studies of stochastic or quantum systems focus on the spectrum of the generator, e.g., the Hamiltonian of the system. The characteristic polynomial is rarely used even though it contains the same information as the spectrum provided that the system is finite-dimensional. A notable exception is the polynomial-based approach to quantum mechanical perturbation theory by Raghunathan [44]. The significant advantage of the characteristic polynomial over the spectrum is that it is always possible to obtain analytical expressions for the former. This particular feature will be fully exploited in the following. In contrast, analytical expressions for the spectrum, i.e., the roots of the characteristic polynomial, cannot be found in general.
We here consider the characteristic polynomial P ξ (x) = det[xI − L ξ ] of the deformed generator L ξ , with I the identity matrix. Generally, P ξ (x) of degree M can be written as The coefficients a µ (ξ) are given by the sum over the principal minors of order (M − µ) of the deformed generator L ξ [45,46] and depend on the variable ξ, except for the coefficient a M −1 . The general expressions for the coefficients a µ (ξ) are unwieldy, but can be readily calculated for system dimensions of interest. Particularly simple exceptions are Note that each generator L ξ has a unique characteristic polynomial P ξ (x), whereas to a given generic polynomial there may correspond none, one or several generators L ξ . In its factored form, P ξ (x) yields the entire spectrum σ ξ of the generator, including λ(ξ), and it thus contains more information than the cumulant generating function G(ξ). Consequently, generators L ξ with identical characteristic polynomials P ξ (x) produce the same counting statistics and will be referred to as being equivalent. The characteristic polynomial is therefore perfectly suited for studying symmetries with respect to transformations of L ξ that leave the counting statistics unchanged. Specifically, P ξ (x) is invariant under arbitrary similarity transformations of the generator L ξ → T −1 L ξ T . However, the structure of L ξ imposed by Eq. (2) is not preserved under similarity transformations, i.e., the transformed matrix is generally not a valid generator. Structure-preserving transformations are, for instance, permutations applied to the matrix elements of L ξ . These transformations preserve the characteristic polynomial P ξ (x) and the structure of L ξ .
It is instructive to consider the example of a stochastic system with unidirectional transitions between nearest-neighbours, i.e., κ i+1,i > 0 for i = 1, . . . , N with periodic boundary conditions and all other rates zero. In this case, the generator L cl ξ is a triangular matrix, except for a single element, and has the characteristic polynomial Clearly, P cl ξ (x) is invariant under any permutation of the rates κ i+1,i and hence there are N! equivalent generators related by permutations which yield identical counting statistics. In addition, the statistics of the observations is independent of the position of the detector.

Reconstruction of characteristic polynomials from cumulants
The characteristic polynomial P ξ (x) plays an important role in ICS since it can be reconstructed from a finite number of cumulants. Reconstructing the full generating function G(ξ) from cumulants, for example, by exploiting Eq. (4), may be possible for a few special cases. To our knowledge there is however no general method available.
To start with, we parametrize P ξ (x) in a more convenient way. In view of the local dependence of the cumulants c ν on the eigenvalue λ(ξ) according to Eq. (4) we constrain our analysis to the properties of P ξ (x) in the vicinity of ξ = 0. We consider the Taylor expansion of the coefficients a µ (ξ) around ξ = 0, which yields with the shorthand notation a Since a µ (ξ) depends on ξ only through the factor e ξ in the generator L ξ we find that a The characteristic polynomial P ξ (x) in the vicinity of ξ = 0 is thus fully determined by the pair of polynomials P ξ (x)| ξ=0 and ∂ ξ P ξ (x)| ξ=0 , or equivalently parametrized by 2(M − 1) of the coefficients a µ and a ′ µ , with µ = 0, . . . , M. Explicitly, the polynomials P ξ (x) ξ=0 and ∂ ξ P ξ (x) ξ=0 are specified by the two sets of Next, the direct relation between the cumulants c ν and the characteristic polynomial P ξ (x) is established. We note that the equality P ξ [λ(ξ)] = 0 holds by definition of the characteristic polynomial. By repeatedly taking the total derivative of this equality with respect to ξ and evaluating it at ξ = 0, we can generate the (infinite) set of equations We observe that Eqs. (9) relate the quantities a µ , a ′ µ and c ν , taking account of the relations c ν = ∂ ν ξ λ(ξ)| ξ=0 for ν ≥ 1 and λ(ξ)| ξ=0 = 0. For instance, the first three functions h ℓ (a µ , a ′ µ , c ν ) are given by As can be seen from Eq. (8) each function h ℓ is linear in the coefficients a µ , a ′ µ . Moreover, considering h ℓ as a function of a µ and a ′ µ , we observe that only h ℓ with ℓ ≥ µ depend on the coefficient a µ , and consequently all h ℓ with ℓ = 1, . . . , M are linearly independent, regardless of the values of the c ν . Lastly, Eqs. (9) are inhomogeneous for ℓ ≥ M because constant terms are generated by the term λ M (ξ).
Using Eqs. (9) we can reconstruct the characteristic polynomial P ξ (x) from the first 2(M −1) cumulants c ν . To this end, we choose the first 2(M −1) equations h ℓ (a µ , a ′ µ , c ν ) = 0 with the cumulants c ν as fixed arguments and solve the resulting linear system for the coefficients a µ , a ′ µ . Aside from exceptional cases, the linear system has a unique solution for the 2(M − 1) coefficients a µ and a ′ µ , thereby yielding a unique P ξ (x). We stress that the reconstruction of the characteristic polynomial P ξ (x) does not guarantee the existence of a generator L ξ compatible with the first 2(M − 1) cumulants c ν . However, if the cumulants c ν indeed result from an evolution governed by a generator L ξ then the reconstructed P ξ (x) is the unique characteristic polynomial of L ξ .
From the reconstructed charateristic polynomial P ξ (x) all cumulants of higher order c ν , with ν > 2(M − 1), can be found. One way of doing so is to evaluate λ(ξ), equivalent to the generating function G(ξ), and subsequently use Eq. (4). We now introduce a direct analytical method, which again is based solely on P ξ (x) and thus avoids the evaluation of G(ξ). Instead of solving Eqs. (9) for the coefficients a µ , a ′ µ with fixed cumulants c ν , we assume that the coefficients are known and solve Eqs. (9) recursively for the cumulants c ν . We illustrate the procedure by evaluating the first three (scaled) cumulants, the average, variance and skewness: where cumulants c ν are expressed in terms of a µ , a ′ µ and cumulants c η with η < ν. In summary, we have two complementary methods at our disposal, which enable us (i) to reconstruct P ξ (x) from a finite number of cumulants and (ii) to find all cumulants from a given P ξ (x). The characteristic polynomial P ξ (x) thus completely replaces the cumulant generating function G(ξ). This remarkable property also suggests an alternative to the standard methods of FCS: Evaluating (zero-frequency) cumulants within FCS usually requires the calculation of the eigenvalue λ(ξ) or the regular part of the resolvant of L ξ [31]. In this traditional way, analytical results can only be obtained for small system dimensions. In contrast, our polynomial-based procedure is efficient and direct: It suffices to calculate P ξ (x) for a given L ξ , which yields the coefficients a µ ≡ a µ (ξ)| ξ=0 and a ′ µ ≡ ∂ ξ a µ (ξ)| ξ=0 , and then to use Eqs. (11). Analytical results for all cumulants c ν can be obtained for any generator L ξ regardless of the system dimension, including L ξ with more general Lindblad operators than in Eq. (2). In particular, the Fano factor is readily found in this way. In general, calculating cumulants recursively seems to be more efficient than finding the full cumulant generating function G(ξ). Recursive schemes for calculating cumulants of high orders have also been exploited in traditional FCS [29][30][31].

Tests for non-classicality, non-Markovianity and system dimension
The reconstruction of the characteristic polynomial makes it possible to design tests for distinguishing between different types of evolutions and system dimensions. The tests are based on the observation that for a finite-dimensional Markovian system, with the Lindblad evolution according to Eq. (2), the number of independent cumulants is finite. Indeed, as shown previously, P ξ (x) can be reconstructed from the first N p = 2(M − 1) cumulants and serves as a cumulant generating function. Thus, the number of independent cumulants is at most N p . The specific value of N p depends on the assumptions about the underlying model Table 1. We design several tests based on the fact that the number of independent cumulants is at most N p . If more than N p independent cumulants are observed in an experiment we can discard the hypothesis under given prior assumptions.

Prior Assumptions Hypothesis
Markovian, Dimension N Classicality Classical, Dimension N Markovianity Quantum, Dimension N Markovianity Classical, Markovian Dimension N Quantum, Markovian Dimension N as regarding classicality, Markovianity and system dimension. If now one of the assumptions is considered on the level of a hypothesis and the number of independent cumulants exceeds the upper limit N p , we can reject this initial hypothesis. As a concrete example, consider a system with two dimensions and Markovian dynamics ('assumptions'). We want to test whether the dynamics is classical ('hypothesis') or not. The system may describe a resonance fluorescence experiment with trapped ions [1] or nitrogenvacancy centers [2], in which the counting statistics of emitted photons is analyzed. Suppose the available cumulants arec 1 ,c 2 andc 3 , where the tilde marks quantities obtained from measurements. Givenc 1 andc 2 , and under the hypothesis that the dynamics is classical we find from Eqs. (10) for the classical case with M = Ñ P cl The predicted classical value of the third cumulant c cl 3 =c 1 + 3c 2 (c 2 /c 1 − 1) then follows from Eq. (11). If the predicted and measured value differ, i.e., c cl 3 =c 3 , then the dynamics of the system is necessarily non-classical. This conclusion relies on the prior assumption that the system is two-dimensional and Markovian, but the argument can be immediately adapted to different combinations of assumptions and hypotheses, as summarized in Table 1. The generic structure of the test is as follows: (i) The first N p + 1 cumulants of the counting statistics are measured. (ii) From the first N p cumulants the characteristic polynomial P ξ (x) is reconstructed in accord with the prior assumptions (classicality, Markovianity or system dimension). (iii) The cumulant c Np+1 is predicted from P ξ (x) and compared to the measured value of c Np+1 .
(iv) If the predicted and measured value of c Np+1 differ then the measured cumulant is independent from the lower-order cumulants. Therefore the dynamics of the system cannot be generated by L ξ and the hypothesis is discarded.
Inevitable experimental uncertainties lead of course to probabilistic rather than sharp test results; however, the result can always be corroborated by measuring and comparing more high-order cumulants and/or using a longer time base for the measurement. Note that the test does not verify the existence of a generator L ξ compatible with the first N p cumulants; this problem is addressed in the next section. As a consequence, the predicted and measured values of c Np+1 can differ for two reasons: Either the generator L ξ compatible with the first N p cumulants makes a false prediction about the measurement result, or the generator does not exist, also leading to a false prediction. In fact, the existence of a characteristic polynomial compatible with measurements already imposes restrictions on the cumulants. We see from Eq. (13), for instance, that for the classical two-state system it is required thatc 1 =c 2 . Moreover, as pointed out in Ref. [24], one can only provide a lower bound on the unknown dimension of the system. In our case, the test reveals that the system dimension must be at least N p +1 and consequently we cannot test for arbitrary dimensions N.

Reconstruction of generators from characteristic polynomials
Even though P ξ (x) contains essential information about the system, it is still desirable to reconstruct the generator L ξ from cumulants. By using the generator together with Eq. (1) we can determine the full time evolution of the system and the steady state ρ s . The reconstruction furthermore allows us to find restrictions imposed on the cumulants that warrant the existence of a generator.
For the reconstruction we have to solve the following inverse problem: Find the values of the parameters entering the structured generator L ξ such that it reproduces the observed cumulantsc ν . The necessary and sufficient condition for L ξ to generate the cumulantsc ν is whereP ξ (x) is the reconstructed characteristic polynomial and P ξ (x) is the characteristic polynomial of L ξ . We recall that while P ξ (x) is invariant under similarity transformations, the structured generator L ξ is not, which restricts the number of solutions to the inverse problem considerably. Let us outline the direct analytical approach to solving the inverse problem. We exploit the fact that analytical expressions are available for P ξ (x) and hence for the coefficients a µ , a ′ µ . The latter are multivariate polynomials of degree M − µ in the parameters of the model; the parameters are collected into the set S for convenience. Equation (14) is thus equivalent to the system of 2(M − 1) polynomial equations From a geometrical point of view, each equation defines an algebraic variety in the parameter space corresponding to S, and the solution of Eqs. (15) is given by the intersection of these varieties [47,48]. It is not necessarily the case that all polynomial equations (15) are algebraically independent, which can be checked, in principle, by using the standard methods of algebraic geometry [47,48]. To take this into account we introduce N i ≤ 2(M − 1), the number of algebraically independent equations. If the number of parameters |S| is larger than N i , the generator L ξ cannot be determined by our approach. Otherwise, if |S| ≤ N i , we can determine all structured L ξ compatible with Eq. (14) in two steps: First, we select |S| independent equations from Eqs. (15) and solve the resulting polynomial system, which reduces the original solution space to a finite discrete set, denoted B. Subsequently, we discard all solutions in B that do not satisfy the full set of Eqs. (15), which finally leaves us with the solutions of the inverse problem.
The first step results in the polynomial system a α (S) =ã α , where the indices α, β ∈ {0, . . . , M − 1} are used to select |S| of the N i independent polynomial equations. In practice, it is advisable to select polynomials with the smallest possible degrees for solving Eqs. (16). At this point, we invoke Bézout's theorem [47,48], stating that Eqs. (16) generally have a finite number of solutions, forming the set B. More precisely, if we admit complex solutions, the number of distinct solutions is at most k d k , where d k is the degree of each polynomial in Eqs. (16). This upper bound scales as |S|! and therefore increases faster than exponentially with the system dimension N. The number of solutions can however be significantly reduced by initially restricting the parameter space, most importantly, to real-valued solutions for stochastic systems.
After selecting from B all solutions that solve the full set of Eqs. (15) we are left with zero, one or a finite number of solutions. No valid solution indicates that the underlying model is not compatible with the observed cumulantsc ν , e.g., false assumptions are made about the dimension N. More than one solution is found if several equivalent generators exist, which by definition cannot be distinguished by the detector. Additional information about the system, such as state occupation probabilities, can be used to single out a unique generator L ξ . Alternatively, several independent detectors can be employed to reduce the number of solutions.

Embedding of the classical into a quantum model
The dimension of L ξ is an important factor in the reconstruction because it sets the number of conditions in Eqs. (15) and hence an upper bound for N i . The number of conditions scales as N for classical systems and N 2 for quantum systems, whereas the number of parameters |S| typically scales as N 2 in both cases. Therefore, ICS applied to classical systems seems to be limited to small dimensions because of the condition |S| ≤ N i .
We can avoid this problem to some extent by treating the classical system as quantum with a trivial Hamiltonian H ≡ 0 in Eq. (2), thereby embedding the classical into a quantum model. The resulting quantum generator L qm ξ , according to Eq. (2), for an actual classical system is sparse and block diagonal, i.e., L qm ξ separates into the block L cl ξ and the diagonal block L coh ξ for the coherences, as exemplified by the generator in Eq. (27). The block L coh ξ contains elements stemming from the terms − 1 2 {L † ij L ij , ρ} in Eq. (2), which in the longtime limit destroy all coherences ρ i =j that may exist initially. As a consequence, the steadystate evolutions resulting from the generators L qm ξ and L cl ξ are the same, and the steady-state coherences are identically zero.
The point of interest for the reconstruction is of course the relation between the classical characteristic polynomial P cl ξ (x) and quantum mechanical polynomial P qm ξ (x) obtained from the embedding. The characteristic polynomial of the block-structured generator L qm ξ factorizes as P qm The additional conditions we gain from treating the classical system as quantum mechanical therefore originate from the polynomial P coh ξ (x), which can be made explicit by writing Eq. (14) in two parts as P cl ξ (x) =P cl ξ (x) and P coh ξ (x) =P coh ξ (x). In practice, starting from 2(N − 1) measured cumulants of the classical system we reconstructP cl ξ (x), and in turn generate the first 2(N 2 − 1) cumulants. Using these 2(N 2 − 1) cumulants we then proceed as for a quantum system to find the rates κ ij . In this way, the embedding strategy extends the use of ICS to classical systems of larger dimension.

Practical reconstruction examples
We illustrate the reconstruction of generators with concrete examples. For the stochastic twostate system, we present the general solutions of Eqs. (16) in terms of the measured cumulants c ν . As further examples we consider a laser-driven atomic system and the Michaelis-Menten kinetics of enzymatic reactions. In the latter cases, we produce the first 2(N 2 − 1) cumulants for a fixed set of parameters and afterwards reconstruct generators L ξ compatible with these cumulants. Incidentally, this application of ICS constitutes a general method for either verifying the uniqueness of a generator or revealing symmetries of the system that are not immediately apparent from the characteristic polynomial P ξ (x). An example for the characteristic polynomial approach in connection with traditional FCS will be presented elsewhere [49].

Stochastic two-state system
We first reconsider the classical two-state system with the detector at the transition |1 → |2 . The corresponding generator reads L cl ξ = −κ 21 κ 12 e ξ κ 21 −κ 12 (18) and the characteristic polynomial is The reconstructed characteristic polynomialP ξ (x) is given by Eq. (13). To find the parameters S = {κ 21 , κ 12 } of the generator L cl ξ we have to solve Eq. (16), i.e., the polynomial system with the solutions {κ 21 , κ 12 } = {κ + , κ − } and {κ 21 , The fact that we obtain two solutions is in agreement with Bézout's theorem and reflected in the symmetry of the characteristic polynomial P cl ξ , which is a special case of Eq. (6). It follows from Eq. (21) that a classical generator L cl ξ exists for cumulants restricted to the regimec 1 >c 2 ≥ 1 2c 1 . In particular, there is no generator L cl ξ corresponding to an observed super-Poissonian (c 2 >c 1 ) counting statistics.

Laser-driven atom with spontaneous decay
We next consider a three atomic states in a Λ-type configuration [50], where the counting statistics is obtained by observing emitted photons (cf. Fig. 2). The three states |1 , |2 and |3 are connected by coherent transitions, parametrized by the Rabi frequencies Ω ij . We assume that the on-site energies of the states are negligible, as is the case for near-resonant laser-driving, so that the Hamiltonian H has the form Spontaneous decay is modelled by the stochastic transitions |3 → |1 and |1 → |2 with decay rates κ 13 and κ 21 , respectively, where the latter is monitored by the detector. The generator of the evolution of the density matrix ρ = {ρ 11 , ρ 12 , . . . , ρ 33 } reads with δ = κ 21 + κ 13 , from which the (albeit large) analytical expression for the characteristic polynomial P qm ξ (x) is readily found.   (23) depend only on the square of the Rabi frequencies, which results in the sign symmetry with eight equivalent generators. The example illustrates that by observing the spontaneous decay of state |1 we can determine the magnitude of all Rabi frequencies Ω ij and spontaneous decay rates κ ij of the atomic system.

Michaelis-Menten kinetics with fluorescent product molecules
Finally, we apply ICS to enzymatic reactions that are described by the Michaelis-Menten kinetics [7][8][9]. The kinetics of a single enzyme can be modelled by a stochastic three-state system, where the states correspond to empty enzyme (E) ≡ |1 , the enzyme-substrate complex (ES) ≡ |2 and the enzyme-product complex (EP) ≡ |3 [51]. The enzymesubstrate binding is a reversible process while the other transitions are assumed to be unidirectional (cf. Fig. 2). The counting statistics is obtained from monitoring the transition (EP) → (E) through the detection of single fluorescent product molecules [7].
The sum of the rates κ ij is identical to the coefficient a N −1 of P cl ξ (x), which holds for any stochastic system and provides a convenient consistency check. We see that in this specific case the counting statistics of the observed product molecules is not sufficient to determine the rates κ ij uniquely. Occupation probabilities would nevertheless provide the possibility to discriminate between the two sets of rates.

Conclusions
We have addressed the problem of finding the properties of a finite-dimensional stochastic or open quantum system from the counting statistics of observable time-discrete events, a scenario encountered in many experimental situations. The first crucial step toward the solution of this newly posed problem was the reconstruction the characteristic polynomial of the deformed generator from a finite number of cumulants by merely solving a linear system. It was moreover shown that the cumulant generating function can be replaced by characteristic polynomial, from which all cumulants can be determined recursively.
By exploiting the fact that only a finite number of the cumulants are independent we have proposed cumulant-based tests for identifying non-classicality, non-Markovianity and a lower bound for the system dimension. In particular, the non-classicality of a Markovian system with dimension N can be identified by measuring the first 2N − 1 cumulants. In contrast to specific state preparations required for similar tests, we only require the system to be in the steady state as we utilize information that naturally leaks out of the system.
As the second important step of ICS, we have suggested a direct analytical approach to the reconstruction of the generator. While perfectly adequate for systems of small dimension this approach requires the solution of potentially large polynomial systems, currently a limiting factor of ICS. Formulated in terms of the spectrum σ ξ , the reconstruction of the generator falls into the class of structured inverse eigenvalue problems (SIEP) [10,52]. In order to apply ICS to larger systems and to cope with inevitable measurement errors we plan to develop efficient and robust numerical methods to solve the polynomial system in Eq. (16) or the equivalent SIEP.
The results of ICS also shed light on the established theory of FCS. It was shown that different generators, i.e., different sets of system parameters, can lead to identical observed counting statistics. Our reconstruction procedure suggests that these equivalent generators form finite discrete sets. The problem of non-uniqueness, an essential aspect of counting statistics, has so far not been directly addressed in the framework of FCS. The characteristic polynomial approach to calculating (zero-frequency) cumulants according to Eqs. (11) offers in addition a potentially useful alternative to the traditional methods of FCS.
Measuring cumulants of high order is experimentally demanding, but in principal a problem of acquiring sufficient data to determine the cumulants with the required precision. Cumulants up to 15th order can be measured in electronic currents through mesoscopic systems [5] and high order cumulants have been recently used to characterize 87 Rb spin ensembles prepared in non-Gaussian states [53]. Interesting stochastic systems to which ICS can be applied are abundant, for example, data transfer in computer networks, traffic problems and biological processes. We believe that ICS as an analytical tool will contribute to a better understanding of these systems.