Quantum Estimation via Sequential Measurements

The problem of estimating a parameter of a quantum system through a series of measurements performed sequentially on a quantum probe is analyzed in the general setting where the underlying statistics is explicitly non-i.i.d. We present a generalization of the central limit theorem in the present context, which under fairly general assumptions shows that as the number $N$ of measurement data increases the probability distribution of functionals of the data (e.g., the average of the data) through which the target parameter is estimated becomes asymptotically normal and independent of the initial state of the probe. At variance with the previous studies [M. Gu\c{t}\u{a}, Phys. Rev. A 83, 062324 (2011); M. van Horssen and M. Gu\c{t}\u{a}, J. Math. Phys. 56, 022109 (2015)] we take a diagrammatic approach, which allows one to compute not only the leading orders in $N$ of the moments of the average of the data but also those of the correlations among subsequent measurement outcomes. In particular our analysis points out that the latter, which are not available in usual i.i.d. data, can be exploited in order to improve the accuracy of the parameter estimation. An explicit application of our scheme is discussed by studying how the temperature of a thermal reservoir can be estimated via sequential measurements on a quantum probe in contact with the reservoir.


I. INTRODUCTION
Seeking the most efficient way to recover the value of a parameter g encoded in the state ρ g of a quantum system is the fundamental problem of a branch of quantum information technologies [1], which goes under the name of quantum metrology [2,3]. It goes without mentioning that this topic has applications in a variety of different research areas, ranging e.g. from the interferometric estimation of the phase shifts induced by gravitational waves [4], high-precision quantum magnetometry [5,6], to remote probing of targets.
In the standard approach one typically focuses on the case where several (say N ) identical copies of ρ g are available to experimentalists, who can hence rely on the statistical inference extracted from independent and identically distributed (i.i.d.) measurement outcomes to estimate the value of g. This scenario is particularly well formulated by those configurations where the unknown parameter g is associated with some black-box transformation Λ g (say a phase shift induced in one arm of an interferometric setup) which acts on the input state ρ 0 of a probing system (say the light beam injected into the interferometer) yielding ρ g = Λ g (ρ 0 ) as the output density matrix to be measured, with such test repeated N times to collect data {s 1 , . . . , s N }. See Fig. 1(a). In this context, the ultimate limits on the attainable precision in the estimation of g, optimized with respect to the general detection strategy, can be computed, resulting in the socalled quantum Cramér-Rao bound, which exhibits the functional dependence upon ρ g via the quantum Fisher information. See e.g. Refs. [2,3,[7][8][9][10][11][12].
In many situations of physical interest, however, the possibility of reinitializing the setup to the same state is not necessarily guaranteed. In the present study we  1. (a) The standard strategy for estimating a parameter g of a quantum system, where measurement data {s1, . . . , sN } are collected by independent and identical experiments. Every time the experiment is performed, the system is reset to some specific known initial state ρ0. (b) The sequential scheme for estimating a parameter g of a quantum system, where the measurements are performed sequentially to collect data {s1, . . . , sN } without resetting the state of the system every after the measurement and the initial state ρ0 can be arbitrary. are going to consider a different scheme, in which a single probing system undergoes multiple applications of Λ g while being monitored during the process without being reinitialized to the same input state. See Fig. 1(b). The data {s 1 , . . . , s N } collected by such sequential measurements will be non-i.i.d. in general. Still, we are able to estimate the target parameter g from the data under certain conditions. We will see that the property of the channel describing the process is important. The idea is to let the probing system forget about the past by the mixing of the channel [13][14][15][16] (the channel being intrinsically mixing or designed to be mixing), which clusters the data and allows the central limit theorem to hold for appropriately chosen functionals of the data. This sequential scheme is suited to account for estimation procedures where one aims to recover g via a sequence of weak measurements that slightly perturb the probe. In particular it can be adapted to study physical setups where the probing system is a proper subset of a many-body quantum system which is directly affected by the black-box generator Λ g (an explicit example of this scenario will be analyzed in the final part of this paper).
Various schemes for quantum parameter estimation based on repetitive or continuous measurements have been studied: see e.g. [17][18][19][20][21][22][23][24][25]. Among them, analogous setups were analyzed in Refs. [19,24], where the problem was formalized in terms of quantum Markov chains. Specifically in Ref. [19] it has been shown that, under rather general assumptions, the statistics of the associated estimation problem converges asymptotically to a normal one, generalizing the similar results which were known to apply to purely classical settings [26][27][28]. In the present paper we first provide an independent derivation of the previous result [19] via a diagrammatic approach to compute the leading-order contributions to the moments of the associated estimating functional of the data {s 1 , . . . , s N }, i.e., the moments of the average S = 1 N N i=1 s i . This approach allows us to prove the central limit theorem including other estimating functionals capturing the correlations among different measurement outcomes, e.g., C = 1 The asymptotic normality of the empirical measure associated to chains of subsequent measurement outcomes is proved in Ref. [24], but in contrast to this previous work we provide explicit formulas which allow us to evaluate the elements of the covariance matrix of the normal distribution of the variables S and C . Moreover, we point out that the inclusion of the correlations C for estimation, which do not contain any useful information in the usual i.i.d. data, can help improve the accuracy of the estimation. This result, while not conclusive, is a preliminary (yet nontrivial) step towards the determination of the ultimate accuracy limit attainable in the non-i.i.d. settings.
This paper is organized as follows. In Sec. II we introduce the notation and recall some basic mathematical facts which will be used in the paper. The non-i.i.d. estimation model is then presented in Sec. III. In Sec. IV we focus on the simplest estimating functional S of the measurement data, and prove its asymptotic normality under the assumption of mixing of the process. The central limit theorem is generalized to include the correlations C and their role in the estimation problem is addressed in Sec. V. An explicit example is then presented in Sec. VI, where we discuss the estimation of the temperature of a thermal reservoir via local measurements on a quantum probe in contact with the reservoir. Conclusions and perspectives are summarized in Sec. VII, while some technical elements are presented in the Appendices.

II. NOTATION AND MATHEMATICAL BACKGROUND
In this section we introduce the notation and recall some basics facts on the theory of quantum channels.
A. Quantum Ergodic/Mixing Channels Quantum channels are completely positive and tracepreserving (CPTP) maps, transforming density operators to density operators [29][30][31]. Every CPTP map E admits at least one fixed point [13,14], namely, a stationary state ρ * , which is hermitian, positive-semidefinite, and of unit trace. In other words, the fixed point ρ * is an eigenstate of the map E belonging to its unit eigenvalue 1.
If the fixed point ρ * is unique, the quantum channel E is called ergodic [14][15][16]. It implies where E n = E • · · · • E denotes n recursive applications of the channel E, and the convergence is in the superoperator norm with corrections whose leading order scales as 1/N . Moreover, if the fixed point ρ * is unique and the unit eigenvalue 1 is the only peripheral eigenvalue (eigenvalue of unit magnitude), the quantum channel E is converging as and is called mixing, with the convergence being as in (2.2) [13][14][15][16]. Mixing implies ergodicity, but the converse is not necessarily true. As commented above, the fixed point ρ * of a quantum channel E is an eigenstate of E belonging to its unit eigenvalue 1. In a matrix representation of E, it is a "right eigenvector." The corresponding "left eigenvector" belonging to the same eigenvalue can be different from the right eigenvector in general. For a quantum channel E, the trace Tr is a left eigenvector belonging to the unit eigenvalue 1, since the quantum channel E is trace-preserving, Tr{E(ρ)} = Tr ρ. (2.4) Let us hence write the fixed point ρ * and the trace Tr in the vectorized notation as ρ * ↔ |ρ * ), Tr ↔ (1|, (2.5) respectively, and a couple of eigenvalue equations for the unit eigenvalue 1 read More explicitly, given any complete set of orthonormal basis states {|n } n of the system, an operator A = n,n A nn |n n | is vectorized by |A) = n,n A nn |n ⊗ |n [32]. The trace (1| = n n| ⊗ n| is (the hermitian conjugate of) the vectorized version of the identity operator: that is why it is denoted by (1|. In addition, the inner product (A|B) = Tr{A † B} is the Hilbert-Schmidt inner product. In this representation, the quantum channel E ↔ m,n,m ,n E mn,m n |m m | ⊗ |n n | is a matrix with the matrix elements E mn,m n = m|E(|m n |)|n in the original representation, and the application of a map E is expressed by the multiplication of the corresponding matrix. By abuse of notation we use the same symbol E for its matrix representation.
The eigenvectors belonging to different eigenvalues are orthogonal to each other and normalized as Note that the matrix E might not be diagonalizable but is cast in the Jordan canonical form in general [13,14]. In this paper, ergodic or mixing channels will play a central role. The unit eigenvalue 1 of such a channel E is not degenerated by definition, and the ergodic/mixing channel E can always be decomposed as where is the eigenprojection belonging to the nondegenerate unit eigenvalue 1 of E, and the remaining part E (which is not CPTP) is built on the eigenvectors {|u n )} n =0 and {(v n |} n =0 belonging to the eigenvalues λ n different from 1. By construction E is orthogonal to P * , i.e., P * E = E P * = 0. Moreover, since E does not admit a unit eigenvalue 1, the inverse (1 − E ) −1 exists, and we have with Q * = 1 − P * , which allows us to prove the convergence to the fixed point ρ * in (2.2). If the channel E is not only ergodic but also mixing, the spectral radius of E is strictly smaller than 1, and we get which proves (2.3). If the channel E is ergodic but not mixing, E admits a peripheral eigenvalue, and E N does not decay: we lose the convergence (2.3), but the averaged channel converges as (2.2). Note again that E might not be diagonalizable if some of the eigenvalues λ n of E are degenerated, but it is not a problem for the convergence: see [13][14][15][16].

B. Measurement and Back-Action
We recall that in quantum mechanics the most general detection scheme can be formalized in terms of positive operator-valued measure (POVM). See e.g. Ref. [29]. Expressed in the superoperator language this accounts to assigning a collection M = {M s } s of trace-decreasing channels M s describing the statistics of the measurement and the back-action on the probed system. In particular, given ρ the density matrix of the system before the measurement, the probability of getting outcome s by the measurement M is given by obtained by summing over all possible values of s, is CPTP and describes the evolution of the system when no record of the measurement outcome is kept. We also notice that given D and E two CPTP maps, the set of channels M s = EM s D defines a new POVM measurement M = {M s } s , where immediately before and after the measurement M one transforms the state of the probed system through the actions of D and E, respectively. Finally we observe that given M = {M r } r and N = {N s } s two POVMs, the operator (N s • M r )(ρ) represents the conditional (not normalized) state obtained when the measurements M and N are performed on a system in the state ρ yielding measurement outcomes r and s, respectively, the associated probability given by p(r, s|ρ) = (1|N s M r |ρ).

III. SEQUENTIAL SCHEME
The problem we study is the following: we wish to recover an unknown parameter g of a quantum system, which is encoded in the state of a quantum probe via the action of a quantum channel Λ g , Here ρ 0 is the input state of the probe, which (possibly) is initialized by us, while Λ g (ρ 0 ) is the associated output state, on which we are allowed to perform measurement in order to learn about g. In a standard i.i.d. approach [10], one is supposed to perform the same experiment several times collecting i.i.d. outcomes {s 1 , . . . , s N }, from which the value of g is to be extrapolated via some suitable data processing. See Fig. 1(a). More precisely, in every experimental run of such an i.i.d. scheme the probe should be initialized in the same input state ρ 0 and the same POVM measurement M = {M s } s should be performed after Λ g has operated on the probe. On the contrary, in the protocol we are going to discuss here, while we keep performing the same measurement M on the probe, the probe is not reset to ρ 0 after each measurement step. Instead, we just repeat the application of Λ g followed by a measurement many times to get a sequence of outcomes {s 1 , . . . , s N }, whose statistics is not necessarily i.i.d. anymore. See Fig. 1(b). In this scenario, following the framework detailed in Sec. II, the state of the probe undergoes a conditional evolution described by the (not necessarily normalized) density matrix defines the probability of the associated measurement event.
It is worth observing that this mathematical setting includes the i.i.d. scenario as a special case, where Λ g is identified with Λ g • P 0 , with P 0 = |ρ 0 )(1| = ρ 0 Tr{ • } being the map resetting the state of the probe into ρ 0 . Indeed, with this choice the probability (3.3) coincides with the one for the case where the measurements M = {M s } s are performed independently on N copies of Λ g (ρ 0 ), i.e., p(s 1 , . . . , s N |ρ 0 ) = (1|E s N |ρ 0 ) · · · (1|E s2 |ρ 0 )(1|E s1 |ρ 0 ) = p(s N |ρ 0 ) · · · p(s 2 |ρ 0 )p(s 1 |ρ 0 ). (3.4) From (2.14) it follows that the map E obtained by summing E s in (3.2) over all s, is CPTP, ensuring the proper normalization of the probability (3.3). As an additional constraint we will require it to be mixing (in some cases, e.g., in Sec. IV A, however, we will weaken this requirement by imposing E to be just ergodic). This is not a strong assumption, as mixing channels actually form an open and dense set. Under this condition we will be able to prove that the parameter g can be estimated from the single sequence of data {s 1 , . . . , s N } collected by the sequential measurements, irrespective of the initial state ρ 0 [19]. The rough idea is that, thanks to the mixing (2.3), repeated applications of the channel force the quantum system to forget its initial state and at the same time decorrelate the data separated beyond the correlation length, which clusters the data and allows us to define self-averaging quantities as estimating functionals, whose fluctuations diminish as N increases, i.e., the central limit theorem holds.
Inferring g from {s1, . . . , s N } A standard estimating functional of the measured data {s 1 , . . . , s N }, through which one tries to infer the value of g, is the average In Ref. [19] it was noted that under the assumption that the average channel E in (3.5) is mixing the central limit theorem holds for S, and for large N the probability distribution P (S) of S asymptotically becomes a Gaussian peaked at a value S * with a shrinking variance σ 2 /N , which are both independent of the input state ρ 0 of the probe, i.e., The explicit expressions for S * and σ will be provided in (4.4) and (4.12), respectively, in Sec. IV below. This ensures that the quantity S evaluated from the single sequence of measurement outcomes is expected, with a high probability, to be very close to its expectation value S * with a vanishingly small variance σ 2 /N for large N . Therefore, by comparing the observed value of S with the formula for the expectation value S * as a function of g, one can infer the parameter g. It is worth stressing once more that in the sequential scheme the measurement data are not independent of each other. Therefore, it is not trivial whether the central limit theorem holds, which is usually based on i.i.d. data set. The mixing, however, is strong enough to kill the correlations between two data if they are sufficiently far away from each other, and clusters the data, allowing the central limit theorem to hold. Thanks to (3.7) the uncertainty in the estimation of g through the quantity S can be evaluated via the Cramér-Rao bound as [3,[7][8][9][10][11][12]33] where F(g) is the Fisher information of the problem given by (3.9) Accordingly, as long as S * exhibits a nontrivial functional dependence upon g, the Fisher information F(g) increases linearly in N , yielding an estimation error (3.8) which diminishes as δg 1/ √ N [in (3.9) we have omitted the contribution from ∂σ/∂g to the Fisher information F(g) since it does not grow with N ].
It may happen however that the quantity S * does not depend upon g. In such a case F(g) nullifies, signaling that it is impossible to recover g through S [a problem which cannot be fixed by properly choosing the input state ρ 0 of the probe, the asymptotic distribution (3.7) being independent of ρ 0 ]. Nonetheless, even in this particular case, the sequence of data {s 1 , . . . , s N }, which is not i.i.d. in general, can still contain some functional dependence upon g, which can be exploited for the estimation of g. In particular, the aim of the present work is to show that the correlations among the measurement data, which are absent in the usual i.i.d. data, can be used for this purpose. It turns out that, under the same mixing assumption on the channel E that leads to the central limit theorem for S in (3.7), the correlations are also self-averaging and become asymptotically normal for large N , enabling one to estimate g through them. See (5.21) and (5.22) in Sec. V below. Even in the case where S * depends upon g, looking also at the correlations help enhance the precision of the estimation of g, which will be demonstrated in Fig. 10 with the example studied in Sec. VI.
We first present an alternative derivation of the results of Ref. [19], i.e., the asymptotic normality of S, on the basis of a diagrammatic approach in Sec. IV. While our approach is more involved than the elegant perturbative approach taken in Ref. [19], it allows us to generalize the scheme to include the correlations among the measurement data in a straightforward manner to enhance the precision of the estimation. We shall indeed prove the asymptotic normality of variables including the correlation functionals in Sec. V.

IV. STATISTICAL BEHAVIOR OF S
This section is devoted to provide an alternative derivation of the results of Ref. [19], which ultimately leads to the asymptotic normality of S in (3.7). We start in Sec. IV A by proving that under the hypothesis that the average channel E in (3.5) is ergodic the quantity S is self-averaging, converging to a fixed value S * independent of the input state ρ 0 . Then in Sec. IV B we introduce the mixing property and show that under this stronger condition the distribution P (S), which rules the statistics of S, becomes asymptotically normal.

A. Law of Large Numbers by Ergodicity
Consider the expectation value of the quantity S with respect to the probability (3.3) governing the statistical distribution of the measurement outcomes, i.e., where E (1) is defined by Assume here that the channel E is ergodic with unique fixed point ρ * : using (2.11), the right-hand side of (4.1) can be written as 3) The first contribution is the value of S N when the input state ρ 0 of the probe coincides with the fixed point ρ * of E. As stressed by the last identity, it also coincides with the expectation value associated with the i.i.d. measurement on ρ * with the POVM M. The second contribution in (4.3) instead is a correction which scales at most as 1/N for any other choice of ρ 0 . Accordingly in the large-N limit we get In a similar way we can compute the variance of S, obtaining where S 2 N is defined similarly to S N in (4.1), and is the variance of s in the stationary state ρ * as in (4.4). Equation (4.6) shows that the variance shrinks as (∆S) 2 N ∼ 1/N , and the fluctuation of S around S N becomes smaller and smaller as we proceed with the measurements. As a result, the probability of finding a singleshot value S close to its expectation value S N becomes very high. Indeed, Chebyshev's inequality bounds the probability of S deviating from S N as for any positive K. In this way, S is self-averaging: each single S is very close to its expectation value with very high probability. In addition, as shown in (4.5), S N becomes independent of the initial state ρ 0 .

B. Beyond the Law of Large Numbers: Central Limit Theorem by Mixing
We have so far assumed the ergodicity of the channel E: this is necessary and sufficient for the convergence S N → S * in (4.5), and for the shrinking variance (∆S) 2 N ∼ 1/N in (4.6). If we further assume that E is mixing, we can say more. For instance, the third contribution to the variance (∆S) 2 N in (4.6) decays as E N /N [note that the sum over j accumulates to O(N )], i.e., faster than 1/N (it is not guaranteed under the ergodicity, since E N does not decay), and the variance (∆S) 2 N asymptotically becomes independent of the initial state ρ 0 . This is because the mixing makes the system forget the initial state ρ 0 as (2.3) without averaging along the time trace.
Most importantly, if E is mixing, one can prove that the probability distribution of S asymptotically becomes normal, converging to the Gaussian distribution (3.7). The asymptotic normality of S under the mixing condition was proved in Ref. [19]. Here we derive the same result by introducing a diagrammatic approach. Specifically, in the following subsection we shall compute the moments of the variable S − S * , showing that for large N they admit the scaling where x denotes the smallest integer not less than x. In particular for even n we shall see that the leading-order term is given by whereẼ (1) andẼ (2) are defined as in (4.7). These results allow us to conclude that the characteristic function for the scaled variable x = √ N (S − S * ) becomes asymptotically normal in the limit N → ∞. Indeed by direct substitution we have (4.13) Accordingly the central limit theorem holds and in the limit N → ∞ the probability distribution P (x) of x converges to a Gaussian peaked at x = 0 with variance σ 2 , i.e., P (x) → e −x 2 /2σ 2 / √ 2πσ 2 , which in the original variable implies (3.7).

Diagrammatic Approach to Evaluate the Moments of S
The expression for the first moment of S − S * follows from (4.3) and is equal to which scales as 1/N as anticipated. Analogously the second moment is readily obtained from (4.6) by noticing that For future reference we find it useful to rederive it: To simplify this we insert the decomposition of the ergodic channel E given in (2.9), namely, we insert P * = |ρ * )(1| or E in place of E. Notice however that Due to this condition, the places in which we can insert P * are limited. The nonvanishing contributions to the second moment hence read 18) and the direct computation of the summations yields (A1) and (4.6).
It is now clear why the second moment (∆S) 2 N in (4.6) as well as the second moment (S − S * ) 2 N scales as 1/N . There are two cases, as we saw in (4.16): (i) two points δs i and δs j coincide (i = j) and we have a single summation 1 N 2 i ; (ii) two points δs i and δs j do not coincide (i = j) and we have double summations 1 N 2 i =j . In any case, once E k with some power k (= j − i − 1 or i − 1 in the above formula for the second moment) is substituted by P * , a summation accumulates as 1 N → O(1), while the contribution from E does not: recall the geometric series in (2.11), where the contribution from E remains O(1/N ). Thus, the substitution rules for estimating the scaling are: Due to (1|Ẽ (1) |ρ * ) = 0 and the coincidence of δs i and δs j , we can insert at most one P * in place of E for the second moment (S − S * ) 2 N : see (4.18). Therefore, the second moment (S − S * ) 2 N is at most O(1/N ), and so is the variance (∆S) 2 N . Note that the second substitution rule in (4.19) is not valid if the channel E is ergodic but not mixing. Indeed, in such a case, the last term in (4.18) yields O(1/N ), as mentioned in the beginning of this subsection. The rule is safe if E is mixing.
We can generalize the above way of estimating the scaling to higher central moments, but a bit more sophisticated rules are required to check the asymptotic normality: we need to care about not only the scalings but also their coefficients. Anyway, the basic strategy to collect the leading-order contributions is to try to insert P * as many times as possible in place of E avoiding (1|Ẽ (1) |ρ * ) = 0. Another important observation is that the insertion of P * = |ρ * )(1| "breaks" the process into pieces. Recognizing these points, we introduce a diagrammatic way of representing the contributions to the moments.
The nth moment is given by Within the summations over {i 1 , . . . , i n } we relabel the n points {δs i1 , . . . , δs in } in chronological order 1 ≤ i 1 ≤ · · · ≤ i N ≤ N and represent them by n dots "•" lined up in chronological order from right to left. See Fig. 2(a). The right most "•" represents the initial state |ρ 0 ), and a trace (1| is supposed to be at the left end. The points can coincide [i = i +1 , as in the case i = j for the second moment: see (4.16)], while between nondegenerate points , which are to be substituted by P * or E i +1 −i −1 Q * , as we did for the second moment [we need Q * to remove i. When E i +1 −i −1 between two points is substituted by E i +1 −i −1 Q * , we connect the two points by a solid line. ii. When E i +1 −i −1 between two points is substituted by P * , we leave the two points disconnected.
iii. In the case where two or more points coincide, we connect the points by dashed lines. (Note that "•" cannot be connected by a dashed line.) See Fig. 2(b). There are two constraints due to (1|Ẽ (1) |ρ * ) = 0: a. The left most two points are surely connected either by a solid line or by a dashed line, since we cannot insert P * between them due to (1|Ẽ (1) |ρ * ) = 0 with the left most trace (1|. b. Each point (except for "•") must be connected with at least one adjacent point either by a solid line or by a dashed line, since we cannot insert P * on both sides of a point due to (1|Ẽ (1) |ρ * ) = 0.
Then, it is easy to draw the diagrams relevant to the leading-order contributions, with the largest possible number of P * inserted. The relevant diagrams for n = 2, 3, 4 are shown in Fig. 3 [see how the two diagrams for n = 2 correspond to the two leading-order terms in (4.18)]. For each diagram contributing to the nth moment: 1. Assign E i +1 −i −1 Q * to each solid line.
2. Insert P * = |ρ * )(1| for each space between disconnected points. · · · to sum the contributions over all possible distances between nondegenerate points respecting the chronological ordering of the points, with an appropriate coefficient counting how many times such a diagram (the specific ordering of the points) appears in the original full range summation 1 N n · · · exploring all possible orderings of the points. The right coefficient reads n!/m 1 ! m 2 ! · · · , where m i are the numbers of coincident points connected by dashed lines in the relevant diagram and the factors m i ! are to disregard the orderings among the coincident points.
It is easily recognized from Fig. 3 that the maximum number of P * we can insert for the nth moment is given by n 2 (where now x denotes the largest integer not greater than x). Therefore, the substitution rules in (4.19) tell us that the nth moment scales as anticipated in (4.10) (the power of N is obtained by n− n 2 = n 2 ). As discussed in the beginning of this subsection, this is the right scaling for the central limit theorem, and only the even moments (n = 2, 4, 6, . . .) are relevant. An important observation is that the leading-order contributions to the even moments are independent of the initial state ρ 0 , since "•" representing the initial state |ρ 0 ) is always disconnected from the first "•". See Fig. 3.
Let us look more carefully at the fourth moment. The leading-order contributions represented by the diagrams in Fig. 3 read This suggests that by the summations each of the leading-order contributions to an even moment (S − S * ) n N acquires a common factor (binomial coefficient) Indeed, each leading-order diagram for an even moment consists of pairs of points (i 2r−1 , i 2r ) (r = 1, . . . , n/2) connected by dashed or solid lines (see Fig. 3), and its evaluation proceeds two points by two points with the help of the following formulas for r = 1, . . . , n/2 (with a convention i n+1 = N + 1): for pairs of coincident points connected by dashed lines [we have a single sum for each coincident pair: see (4.21)], while for pairs of nondegenerate points connected by solid lines [The actual ranges of the summations in the leadingorder contributions to the even moments are slightly different from those in (4.22) and (4.23), but the corrections are finite and become negligible in the asymptotic regime N n.] In this way, each leading-order diagram for an even moment acquires the binomial coefficient N !/(n/2)!(N − n/2)! with E i +1 −i −1 being transformed into (1 − E ) −1 . This leads us to the following recipe for obtaining the expressions for the leading-order contributions to any even moment directly from the relevant diagrams: 1 . Assign (1|Ẽ (2) |ρ * )/2 to each pair of points connected by a dashed line.
Then, the leading-order contributions to the even moments factorize as shown in Fig. 4 and yield (4.11).

V. USE OF CORRELATIONS
An important difference from the standard strategy for parameter estimation, where independent identical experiments are performed to collect data, is that in the present sequential scheme the correlations among the measurement data are available for estimation. Combining the information attainable from the correlations with that from the average S, the precision of the estimation can be enhanced. The primary motivation of the present paper is to explore this possibility.
For instance, one can compute from a single sequence of N measurement outcomes {s 1 , . . . , s N }, which captures the correlation between two data separated by a distance . In the presence of the correlations among the data, C may depend on the target parameter g in a way that cannot be deduced solely from S. This might provide additional knowledge on how the parameter g is encoded in the process and can enhance the precision of the estimation of g.
In principle, ranges = 1, . . . , N − 1, but recall that the correlation between two data are expected to decay exponentially as increases under a mixing channel: C with greater than the correlation length would not contain useful information. In addition, N should be much greater than so that the number N − of data used to evaluate C is large enough. Therefore, we will require N L ≥ , with L being the maximum we take to estimate the parameter g.
The correlations C are also self-averaging quantities. Moreover, we are able to prove that the central limit theorem holds for the set of quantities X = (S, C 1 , . . . , C L ). First, the expectation value of C is evaluated as which approaches under the ergodicity of the channel E. Then, let us look at the nth moment (5.4) It is an nth-order polynomial of k = (k 0 , . . . , k L ), and is a collection of all the nth moments among X = (S, C 1 , . . . , C L ) as its coefficients. Since we are interested in the asymptotic limit N → ∞, we collect the leadingorder contributions to µ n (k) for large N . The idea to do that is basically the same as that for S: we try to insert P * as many as possible in place of E between points δs i from S − S * and pairs of points

5)
, while if it does not we leave it disconnected from or connect it with its adjacent dot by a solid line depending on whether P * is inserted between them or not. The ranges of the summations · · · exploring all possible distances between dots "•" should be carefully arranged depending on whether the dots represent points δs i or pairs of points δ(s i s i+ ) , and some of the prefactors in 1/N n are replaced by 1/(N − ), but such details become irrelevant in the asymptotic regime N n, L. Then, the analysis goes in the same way as before, the leading-order diagrams are again given by Fig. 3, and the nth moment µ n (k) for an even n asymptotically factorizes pairwise as Fig. 4, where the pair of dots "•" connected by a dashed line or a solid line represents the collection of all pairwise combinations among δs i and δ(s i s i+ ) ( = 1, . . . , L), with some care on the coefficients to distinguish different orderings of the pieces, i.e., give a coefficient 1/2 to the pair connected by the solid line in Fig. 4 and collect contributions with different orderings of the pieces [see the second and third terms in Σ 00 , Σ 0 , and Σ in (5.9)-(5.11) below]. We get withẼ (m) andẼ (1) defined in (4.7) and (5.7), respectively, 14)  (5.19) corresponding to the diagrams in Fig. 5. We provide the complete expressions for the covariances among S and C valid for any (even small) N in Appendix A, whose asymptotic forms coincide with the covariances (5.9)-(5.11) divided by N . This result shows that the set of scaled variables √ N (S − S * ) and √ N (C − C * ) asymptotically become normal in the limit N → ∞. The characteristic function reads where Σ is the (L+1)×(L+1) matrix with its matrix elements given by the covariances in (5.9)-(5.11). The central limit theorem holds, and the probability distribution P (X) of X = (S, C 1 , . . . , C L ) becomes asymptotically Gaussian, peaked at X = X * with a shrinking covariance Σ/N . This ensures that the single-shot values X computed from a single sequence of measurement data well represent their expectation values X * , through which we can estimate a parameter g. The uncertainty δg in the estimation of g is given by (3.8) with the Fisher information which increases linearly in N , and the uncertainty δg diminishes as δg 1/ √ N [in (5.22) we have omitted the contribution from ∂Σ/∂g to the Fisher information F L (g) since it does not grow with N ]. Moreover, this Fisher information F L (g) for the estimation of g through a set of quantities X = (S, C 1 , . . . , C L ) can be greater than the Fisher information F(g) = F 0 (g) given in (3.9) for the estimation of the same g but solely through S. The precision of the estimation can be enhanced by looking at the correlation data C in addition to the average S.
Here we have considered the two-point correlations C as well as the average S. If we incorporate higher-order correlations with more points, the precision of the estimation can be further improved. On the other hand, correlations with too many points would not be helpful, since the number of data used to evaluate such correlations is reduced, and some of the points involved in the correlations are separated beyond the correlation length of the mixing channel supplying no more information than lower-order correlations. It is currently not clear to what extent we can improve the precision of the estimation by looking at higher-order correlations.

VI. EXAMPLE: ESTIMATION OF THE TEMPERATURE OF A RESERVOIR
In this section we analyze an explicit example, where the correlations among the data collected by the sequential measurements would be useful for improving the estimation of a parameter. The setting we consider is related to quantum thermometry, which aims to use lowdimensional quantum systems (say qubits) as temperature probes to minimize the undesired disturbance on the sample (see e.g. Refs. [34,35] and references therein). Specifically we focus on the paradigmatic example with a qubit probe in contact with a thermal reservoir at a finite temperature T . Our goal is to estimate the temperature T of the reservoir by monitoring the relaxation dynamics induced on the qubit, which effectively plays the role of a local "thermometer" (Fig. 6). In our approach we describe the probe-reservoir coupling in terms of the resulting Markovian master equation [30,[36][37][38][39] probe qubit reservoir FIG. 6. We estimate the temperature of a thermal reservoir though measurements performed on a probe qubit in contact with the reservoir. operating on the probe, i.e., where ρ(t) represents the state of the qubit, Ω is the energy gap between the excited |↑ and ground |↓ states of the qubit, and The two relaxation constants γ + (for decay) and γ − (for excitation) are related to the temperature of the reservoir T , respecting the detailed balance condition. For a bosonic thermal reservoir, they are given by [36][37][38][39] γ + = (1 + n th )γ, γ − = n th γ, n th = 1 e Ω/k B T − 1 , (6.3) with k B being the Boltzmann constant. We assume that the parameters Ω and γ (i.e., the characteristics of the thermometer) are known. Estimating the temperature T is then equivalent to estimating while γ = γ + − γ − is a known constant independent of the temperature T . The higher is the temperature, the larger is the decay rate γ β .

A. Standard Strategy
The information about the temperature T , namely, the parameter γ β , is imprinted in the state of the qubit through the dynamics under the influence of the thermal reservoir, i.e., by the action of the quantum channel Λ t which is the solution to the master equation (6.1). Then, the standard strategy to estimate the parameter γ β is (i) to prepare the qubit in a specific initial state ρ 0 , (ii) to let the qubit evolve ρ(τ ) = Λ τ (ρ 0 ) for a certain time τ in contact with the thermal reservoir, and (iii) to measure a specific observable in the state ρ(τ ).
We repeat this experiment N times to collect measurement results, from which we estimate the parameter γ β . For instance, we prepare the qubit in a specific initial state ρ 0 , say in the excited state |↑ , and after a fixed waiting time τ we measure the qubit to check whether it is in the excited state |↑ or in the ground state |↓ . We repeat this process N times, and we estimate γ β from the survival probability of the initial state |↑ after time τ . Our measurement however can be weak and unsharp: here we consider the measurement which provokes the following back-action on the qubit, depending on the outcome of the measurement s. This measurement process can be simulated with a cnot gate [40,41]. The parameter η controls the precision and the strength of the measurement: η = 0 provides the projective measurement, while with η = π/4 the measurement gives totally random results with no disturbance on the measured system. The probability of obtaining the measurement outcome s in the state ρ(τ ) is then given by sin 2 η |↑ ↑| + cos 2 η |↓ ↓| (s = −1) (6.8) are the POVM elements of this measurement. The uncertainty in the estimation is then bounded by the Cramér-Rao inequality [3,[7][8][9][10][11][12]33] δγ β 1 N F (γ β ) (6.9) with the Fisher information given by For the present model, the Bloch vector of the qubit evolves as 11) where σ x = σ + + σ − and σ y = −i(σ + − σ − ). The equilibrium state ρ eq is characterized by σ x eq = σ y eq = 0, σ z eq = − γ γ β , (6.12) namely, The probability distribution of the outcomes of the measurement (6.7) at time τ reads p τ (±1|ρ 0 ) = 1 2 1 ± σ z τ cos 2η , (6.14) and the Fisher information F (γ β ) in (6.10) is estimated to be A larger Fisher information would be attainable by measuring a different observable. The maximum Fisher information one can reach with the optimal measurement is given by the quantum Fisher information [2,3,[7][8][9][10][11][12], with the symmetric logarithmic derivative where V is a 3 × 3 matrix whose matrix elements are given by Notice here that both the Fisher information F (γ β ) in (6.15) and the quantum Fisher information F Q (γ β ) in (6.16) depend on the choice of the initial state ρ 0 . Because of the convexity of the quantum Fisher information, the maximum of the quantum Fisher information (the best estimation) is always achieved by choosing a pure input state ρ 0 = |ψ 0 ψ 0 | [3,42]. Moreover, for the present problem, the ground state of the qubit |ψ 0 = |↓ is the optimal choice, in the sense that the maximum of F Q (γ β ) for a given temperature is achieved with |ψ 0 = |↓ : see Fig. 7, and the temporal behavior of F Q (γ β ) for |ψ 0 = |↓ is plotted in Fig. 8. For this specific initial state ρ 0 = |↓ ↓|, the Fisher information F (γ β ) in (6.15) with η = 0 coincides with the quantum Fisher information F Q (γ β ) in (6.16), for any time τ and   9. The temporal behavior of the Fisher information F (γ β ) given in (6.15) for ρ0 = |↓ ↓| and γ β /γ = 1.5 with different strengths of the measurement η. In the case of projective measurement η = 0, the Fisher information F (γ β ) coincides with the quantum Fisher information FQ(γ β ) given in (6.16) and plotted in Fig. 8.
for any γ β : the projective measurement to discriminate |↑ and |↓ is the optimal measurement. For nonvanishing η > 0 the Fisher information F (γ β ) is reduced, and the weaker is the measurement, the smaller is the Fisher information F (γ β ), as shown in Fig. 9.

B. Sequential Scheme
Let us now turn our attention to the sequential scheme. First, it is important to check whether the channel E defined in (3.5) with E s in (3.2) is mixing. For the present model, the spectrum of E is given by {1, e −γ β τ , e −(γ β /2±iΩ)τ }, (6.19) and therefore, E is mixing for any τ > 0 with a unique fixed point (the eigenstate belonging to the eigenvalue 1) 20) which coincides with the equilibrium state ρ eq in (6.13) of the free relaxation process. This mixing is apparently a direct consequence of the irreversibility of the relaxation process Λ t of the probe qubit. Since E is mixing, the sequential scheme works for the present problem. Let us take the average of the outcomes of a sequence of N measurements, S defined in (3.6), as a quantity through which we estimate γ β . For the present model, its expectation value is computed to be (N ≥ 1) (6.21) and the variance to be As N increases, both become independent of the initial state ρ 0 , and the variance (∆S) 2 N shrinks as 1/N , In other words, S evaluated from a single sequence of measurements almost certainly exhibits a value very close to its expectation value S N , which is a function of γ β . Therefore, by comparing S (obtained via a single experimental run) with its expectation value S N [given by the formula (6.23)], the parameter γ β is estimated with the uncertainty regulated by the variance (∆S) 2 N in (6.24), i.e., with the precision given by the Fisher information F(γ β ) = F 0 (γ β ) in (3.9), , (6.25) which is to be compared with the Fisher information N F (γ β ) with (6.15) by the standard strategy (see Fig.  10 below).
As stressed above, the correlations among the acquired data are also available for the estimation in the sequential scheme. For instance, the two-point correlations C defined in (5.1) can be used to estimate γ β . Their expectation values (for a generic initial state ρ 0 ) are given by ( ≥ 0, N ≥ + 1), (6.26) and their covariances (in the stationary state ρ 0 = ρ * ) by [the complete expression for the covariance C C N − C N C N valid for any , ≥ 1 and N ≥ max( , )+1 (but for ρ 0 = ρ * ) is given in Appendix B]. All the covariances scale as 1/N , and the Fisher information (5.22) increases linearly in N . This ensures that, by comparing the set of quantities (S, C 1 , . . . , C L ) evaluated from a single sequence of measurement data with the set of their expectation values ( S N , C 1 N , . . . , C L N ), one can estimate γ β with the precision given by the Fisher information F L (γ β ) computed by the formula (5.22), which increases linearly in N . It is reasonable to expect that the estimation with the multiple quantities (S, C 1 , . . . , C L ) is better in precision than the estimation solely through the average S, namely, the Fisher information F L (γ β ) (L > 0) is larger than the Fisher information F 0 (γ β ), and the more correlations are incorporated (the larger is the number L), the larger is the Fisher information F L (γ β ). Let us look at two different regimes.

Projective Measurement η = 0
The Fisher informations F L (γ β ) (L = 0, 1, 2) by the sequential scheme are plotted in the five panels in Fig. 10(a) and are compared with the Fisher information F (γ β ) by the standard strategy with ρ 0 = |↓ ↓|, for the case of projective measurements η = 0. In this case, the Fisher information F (γ β ) coincides with the quantum Fisher information F Q (γ β ) in the standard strategy.
Compare first F (γ β ) and F 0 (γ β )/N (per measurement). We observe that the standard strategy provides better estimation than the sequential scheme. Recall here that the input state ρ 0 = |↓ ↓| for F (γ β ) is the optimal for the standard strategy. On the other hand, in the se-quential scheme, the state of the qubit is projected into |↑ or |↓ depending on the outcome of the projective measurement. If it is projected into |↑ by a measurement, it restarts to evolve from this non-optimal state for the next measurement. Not all the steps in the sequential measurements are optimal for the estimation. That is why the sequential scheme cannot beat the standard strategy, in the case of projective measurement.
One can improve the performance of the sequential scheme, by incorporating C 1 for the estimation. Indeed, as is clear from Fig. 10(a), the Fisher information F 1 (γ β ) for the estimation through (S, C 1 ) is greater than the Fisher information F 0 (γ β ) solely through S. Note that no additional resources or experiments are required to incorporate C 1 : one simply needs to carry out additional data analysis to compute C 1 from the data used to evaluate S. In Fig. 10(b), the gain in the Fisher information by incorporating C 1 is shown for different temperatures.
On the other hand, incorporating more correlation data, i.e., C with > 1, does not help improve the estimation. See Fig. 10(a) again. This is because every time one performs measurement the system is reset to a pure state by the projective measurement: there is no correlation between the measurement results separated over two steps. The system simply repeats the same dynamics, jumping between pure states |↑ and |↓ , and the measurement after multiple steps gains no more information than that attainable by the measurement after a single step.

Weak Measurement η > 0
Let us next look at the cases with weak measurements η > 0. As is clear from Fig. 10(c), the sequential scheme       can be better than the standard strategy. In particular, at low temperatures, the Fisher informations F L (γ β )/N (L ≥ 0) by the sequential scheme exceed the Fisher information F (γ β ) by the standard strategy. The reason is the following. In the standard strategy, the weak measurement is performed only once, and the system is reset to the specific initial state ρ 0 for the next measurement. The single weak measurement can acquire less information than a projective measurement, but if it is repeated many times, as in the sequential scheme, the information is accumulated, and better information is gained in our hands. At the same time, the system is gradually projected to one of the eigenstates of the measured observable by the repeated weak measurements [43]. In other words, the repetition of the weak measurements mimics a stronger measurement (closer to a projective measurement). That is why the sequential scheme can be better than the standard strategy, in the case of weak measurement.
It is also clear from Fig. 10(c) that the precision of the estimation is improved by incorporating the correlation data C . The gain in the Fisher information [F L (γ β ) − F L−1 (γ β )]/F 0 by adding a correlation C L to (S, C 1 , . . . , C L−1 ) is shown in Fig. 10(d). The enhancement is reminiscent when the time interval between measurements τ is short, i.e., γ β τ 1. Moreover, the gain exhibits a peak at a smaller τ for a larger L. This is because the two points of each two-point correlation C , separated by steps, should be within the correlation time τ c ∼ 2/γ β [which is ruled by the second largest eigenvalues e −(γ β /2±iΩ)τ of the mixing channel (6.19)], in order for the correlation C to bear useful information.
It appears that the sequential scheme can beat the standard strategy only at low temperatures (small γ β ), but it should be noted that the standard strategy in Fig.  10 assumes the optimal initial state ρ 0 = |↓ ↓|, while in the sequential scheme the system is around the stationary state ρ * of the mixing channel E, which is the thermal equilibrium state ρ eq [see (6.20)]. It would be more appropriate to compare the Fisher informations F L (γ β )/N by the sequential scheme with the Fisher information F (γ β ) by the standard strategy in the large τ limit (which gives the Fisher information with the initial thermal state ρ 0 = ρ eq ).

VII. CONCLUSIONS
The estimation of a parameter encoded in a quantum probe, through a series of measurements performed sequentially on the probe, has been analyzed in a general non-i.i.d. setting. On the basis of a diagrammatic approach we have discussed the conditions under which the central limit theorem holds as the number of mea-surements increases, reproducing the previous results [19] and generalizing them to the case where the correlations among the measurement data are also taken into account in the estimation strategy. Our analysis explicitly shows that the latter strategy can yield a significant advantage over the standard procedure where only the average of the acquired data is considered.
At present however it is not clear whether this is the best strategy one can do: it is indeed possible that different data processing (including the evaluation of higherorder correlations commented at the end of Sec. V) can improve further the attainable accuracy. In the example studied in Sec. VI, the sequential scheme surpassed the standard i.i.d. procedure when we are able to perform only weak measurements, but could not beat the standard procedure when we are allowed to perform strong measurements. A better strategy for the sequential scheme could beat the ultimate precision achievable by the standard strategy. The optimal strategy would require different measurements step by step, or moreover would require quantum-correlated measurements over different measurement probings. The use of entanglement is also an interesting possibility [44]. It is yet to be clarified what is the ultimate accuracy attainable in the sequential scheme for parameter estimation [45].
Recently, quantum metrology in the presence of noise is under intense study [44,[46][47][48]. The mixing property required for the sequential scheme is relevant to noisy channels, and connections with such issue would be interesting to be explored. and whereẼ (m) and (∆s) 2 * are defined in (4.7) and (4.8), respectively,Ẽ (m) are in (5.7) and (5.12), and the other components are given in (5.14)-(5.19).
Appendix B: Covariances among C for the Model In (6.27) in Sec. VI we showed the asymptotic expression for the covariance between C and C for large N for the model. Here we provide its complete expression valid for any (even small) N . In the stationary state ρ 0 = ρ * , the covariances between C and C ( ≥ ≥ 1) are given for N ≥ + by while for + > N ≥ + 1 by (B2)