Continuously monitored quantum systems beyond Lindblad dynamics

The dynamics of a quantum system, undergoing unitary evolution and continuous monitoring, can be described in term of quantum trajectories. Although the averaged state fully characterizes expectation values, the entire ensemble of stochastic trajectories goes beyond simple linear observables, keeping a more attentive description of the entire dynamics. Here we go beyond the Lindblad dynamics and study the probability distribution of the expectation value of a given observable over the possible quantum trajectories. The measurements are applied to the entire system, having the effect of projecting the system into a product state. We develop an analytical tool to evaluate this probability distribution at any time t. We illustrate our approach by analyzing two paradigmatic examples: a single qubit subjected to magnetization measurements, and a free hopping particle subjected to position measurements.


Introduction
Quantum mechanics is a fundamental milestone of the human comprehension of natural world [1].One of its most enigmatic and controversial features is the role played by measurements [1][2][3][4][5].While in a classical perspective measurements are trivially a way of extracting information from a system, in quantum mechanics the meaning of measurements is much more profound [6].When a quantum system is measured, its wave function undergoes a non-deterministic collapse and the system is projected into a specific state [7].Recent technological advancements have enabled increasingly accurate and fast measurements on quantum systems [8,9].These advances have opened up new avenues for exploring the fundamental principles of quantum mechanics [10][11][12] and have led to the development of novel applications in fields such as quantum information processing [13,14] and quantum thermodynamics [15][16][17][18][19][20][21][22][23].However, the dynamics of continuously monitored quantum systems are often difficult to analyze, and there is a need for theoretical approaches that can provide insights into the behavior of these systems Great interest has been shown in measurement-induced criticality, which has emerged as a prominent topic in the study of continuously monitored quantum systems .Indeed, continuous measurements of a quantum system can create a feedback loop that leads to a non-equilibrium steady state.In this steady state, the system can exhibit different phases depending on the strength and type of measurement.For example, in the quantum Zeno phase, frequent measurements can suppress quantum fluctuations, leading to a phase where the system behaves as if it were frozen [50,51].Conversely, when measurements are less frequent or weaker, the system can enter a volume law phase, where the entanglement entropy grows linearly with the system size.Nonetheless, in experiments on measurement-induced criticality, post-selection is often required to reveal the underlying physics hidden in quantum trajectories [52,53].
In this article, we consider quantum systems that are coupled to a measuring apparatus, examining their evolution according to their unitary dynamics, which is interrupted by projective measurements with a fixed rate γ.The measurements are applied to the entire Hilbert space, thus projecting the system onto an eigenstate of the measured observable.To illustrate our approach, we derive analytical results for two paradigmatic examples: a single qubit measuring its magnetization and a free hopping particle measuring its position.We therefore provide an exact method for computing the probability distribution of the expectation value of observables averaged over the set of quantum trajectories.
The approach we present here goes beyond a Lindblad master equation description.Indeed, the Lindblad equation describes the dynamics of a quantum system subjected to quantum jumps by considering the density matrix averaged over all possible outcomes, which is a non-selective state [6,[54][55][56].In the case of a continuously measured system, the Lindblad approach neglects the outcome of the measurements and averages over all possible outcomes.In contrast, our approach is selective, as we keep the information contained in individual quantum trajectories.This allows us to uncover the physics that is neglected in the average state.

Protocol
Let us consider an N -level quantum system described by a time-independent hamiltonian H, whose unitary evolution is governed by U (t) = e −iHt .In addition, all along the evolution, we couple the system to a measuring apparatus that project, with a fixed measurement rate γ, the evolved state in to an eigenstate of the observable such that [A, H] = 0.At time t = 0 the system has been prepared in a fixed state |ψ(0) corresponding to an eigenstate |a 0 of A. We are interested in evaluating the probability distribution of the expectation value of a generic observable O = a o a |a a| commuting with A, i.e.
with t > 0. Here, the over-line is indicating the average over the quantum trajectories, which are labelled by the integer ξ.A sketch of our dynamical protocol is shown in Figure 1.We can split this average by putting apart the contributions for different n clicks of the measuring apparatus.For each n the system is projected in one of the eigenstates of A at times {s j }, with j = 1, . . ., n such that 0 < s where P (0) , P (n>0) are respectively the no-click (n = 0) and click (n > 0) contributions to the probability distribution of ψ ξ (t)|O|ψ ξ (t) .Now we can write down these object as follows and p(x, t|s n , a n ; ...; s 1 , a 1 ; 0, a 0 ) is the conditional probability density of getting the average value x at the final time t, given that the system was found in the eigenstates {|a j } of A at times {s j }.T denotes the time-ordered product.The term γ n e −γt is the Poisson weight associated to the n−measurements events (the factor 1/n! is removed since inside the time-ordered integral the events labelled by 1, 2... n occurs at fixed times s 1 < s 2 < ... < s n ).Since the measurement process is Markovian, we have that p(x, t|s n , a n ; . . .; where the delta contribution can be rewritten as and gives the contribution to the probability due to the last time-interval (s n , t).In addition, the transition probabilities between eigenstates of A read Let us assume that the transition matrix T can be diagonalised, i.e.T(s − s) = V D(s −s) V † with eigenvalues D α,β (t) = d α (t)δ α,β , and time-independent unitary matrix V whose columns are the orthonormal eigenvectors.Notice that, since T is a unitarystochastic matrix, d 1 = 1 is the largest eigenvalue, meanwhile as a consequence of the Figure 1: Scheme of the dynamical protocol.We consider a system prepared in a fixed initial state |ψ(0) = |a 0 , undergoing random projective measurements of an observable A (see Eq. ( 1)).Measurements occur with a fixed rate γ at times s 1 < s 2 < ..., whereas in the intervals (s j , s j+1 ) the system follows an unitary time evolution.Different quantum trajectories are labelled by the integer ξ.The ensemble of states |ψ ξ (t) defines the probability distribution function of ψ ξ (t)|O|ψ ξ (t) (right).
Perron-Frobenius theorem |d α | ≤ 1 for α = 2, . . ., N , where the equality holds in case of degeneracy [57,58].In addition, since a T a,b = b T a,b = 1, the eigenvector associated to largest eigenvalue is simply given by V a,1 = 1/ √ N .Now let us focus on the nontrivial P (n>0) O (x; t).For each order n and fixed states {a 0 , a n }, the sum over all the the possible intermediate outcomes can be rewritten as with starting time s 0 = 0. We thus get where we defined the time-ordered integrals which satisfy the following recursive relations In Eq. ( 10) the sum over n can be safely taken, since s n and a n are dummies variables.Thus, we can introduce I α (s) = ∞ n=1 I n,α (s) which satisfy the following linear Volterra integral equation of the second kind Exploiting the Laplace transform of the integral kernel, i.e.L [d α ](z) = ∞ 0 e −zt d α (t)dt, we can easily solve the previous equation obtaining where L −1 [. . .] is the inverse Laplace transform.As expected, it easy to show that thus identifying We can finally collect all those results, further simplify and finding the following general expression for the probability distribution Let us finally mention that using this expression we can easily compute the k-th moment of the distribution, defined as The first moment (k = 1) is particularly simple, since it does reduce to the expectation value of O over the averaged density matrix |ψ ξ (t) ψ ξ (t)| whose dynamics is fully described by the Lindblad equation (see Appendix C).In particular, using that b T a,b (t − s)V b,α = d α (t − s)V a,α , we easily get which can be further simplified exploiting Eq. ( 13), finally obtaining Let us stress that the simplified expression in Eq. ( 18) applies only for averages of diagonal observables.Every time we are interested in higher moments or probability distribution of non-diagonal operators, the computation have to be carried out from scratch, basically starting from Eq. ( 16).Finally, due to the properties of the transfer matrix T, the stationary limit t → ∞ of the previous average can be easily taken; in fact, only α = 1, with I 1 (t) = e γt , will contribute to the sum over α, leading to o a , where we used the fact that V a,1 = 1/ √ N .As expected from the Lindblad dynamics, it does correspond to the expectation value of the operator O over the infinite temperature state 1 N N a=1 |a a|.

Two-level system
We start by considering a continuously monitored two-level system.Despite their simplicity, two-level systems are the fundamental building block of the most studied quantum many-body systems and therefore understanding their behavior is crucial.The system consists of a single spin-1/2, whose unitary evolution is governed by the following hamiltonian where σ α are the Pauli matrices with α = x, y, z and σ α , σ β = 2i αβγ σ γ .The system is continuously monitored along z with a rate γ (therefore A = σ z ).We will denote with |σ = |+1 , |−1 the eigenstates of σ z with eigenvalue σ = ±1 respectively.We consider a two-level spin which starts from |+1 and we want to evaluate i.e. the probability distribution of the expectation value of O = σ z .It is easy to find that I+sin(Jt) 2 σ x , and V = (σ z +σ x )/ √ 2 whose columns are the eigenvectors of σ x .As expected d 1 = 1 and d 2 (t) = cos(Ωt), where Ω = 2J is the splitting between the energy eigenvalues of the system.Notice that, from these eigenvalues we can derive two recursion relations as in the integral equations (12), then, since d 1 = 1 we immediately get I 1 (s) = e γs .In the following, we will simplify the notation I 2 (s) = I(s).Finally, we notice that D(x, σ n , t − s n ) = δ x − σ n cos(Ω(t − s n )) , from which we can compute the no-click contribution P (0) In order to find the full probability distribution, we can use equation ( 16) to write For a better understanding of the solution, we do not apply the Laplace method to solve the integral equation ( 13) for I(s).Instead, we differentiate it twice with respect to s, which yields the following differential equation for I(s) with the initial conditions I(0) = 1 and İ(0) = γ.This is the equation of motion of an "anti-damped" harmonic oscillator (since γ > 0).The solution can be found easily using standard methods [59], obtaining where we defined the following characteristic frequency Γ = γ 2 − 4Ω 2 .As expected, we find three different regimes depending on the sign of γ 2 − 4Ω 2 .Let us analyze the behavior of the function e −γs/2 I(s).When γ < 2Ω, we have that Γ is an imaginary quantity.Therefore, we find a regime in which the function oscillates with a frequency of Ω 1 − γ 2 /(4Ω 2 ), which is lower than the natural frequency Ω.When γ = 2Ω, we get the critical regime for which the function e −γs/2 I(s) becomes linear, (1 + γs/2).Finally, if γ > 2Ω, it grows exponentially with a rate γ/2 1 − 4Ω 2 /γ 2 smaller than γ/2.Let us finally consider the Dirac-delta function contribution.Solving for x − σ cos(Ω(t − s)) = 0, we easily get Using the properties of the Dirac-delta distribution, we can rewrite it as Since in Eq. ( 22), the integral is taken over a finite interval [0, t], we need to constraint the solutions sk in Eq. ( 25) only to those falling in that region.However, we may relax that condition by explicitly introducing the Heaviside step function θ(x), such that Eq. ( 22) finally reads where sk implicitly depends on x and t as well.In Figure 2, we plot the entire distribution P σ z (x; t).As expected, at early time the probability is highly asymmetric, having support only in the vicinity of x = 1.The no-click term is in fact localised at x = cos(Ωt) but its weight is exponentially suppressed in time.The click contribution is also showing a nontrivial asymmetric evolution.Its behavior strongly depends on γ indeed: in the oscillatory regime (γ < 2Ω), the probability distribution function bounces back and forth between the two extremes x = ±1 of its domain; after many oscillations, the number depending on the value of γ, it is expected to relax toward a symmetric distribution, the typical relaxation time being τ = 2/γ.For γ > 2Ω the probability is not oscillating anymore and the relaxation time becomes τ = 2/(γ −Γ).Interestingly, as the measurement rate is getting higher, the time needs for the magnetisation statistics to reach the equilibrium becomes larger and larger, diverging as τ ∼ γ/Ω 2 .The first moment of P σ z (x; t), namely the magnetisation average x P σ z = e −γt I(t), does coincide with the expectation value over the averaged state (see Appendix B).Nevertheless, our approach allows to easily compute the fluctuations of the magnetization along the trajectories.Indeed, by means of Eqs. ( 21) and ( 22), we can easily get the second moment of the distribution which simplifies to for a two-level system (Ω = 1).Notice that the critical value of γ = 2 leads to the fastest convergence to the equilibrium stationary value of the average magnetization.and which cannot be computed using the a Lindblad approach since it corresponds to ψ ξ (t)|σ z |ψ ξ (t) 2 .In Figure 3, we represents x P σ z , x 2 P σ z as a function of time for different values of the measurement rate γ.
Finally, let us mention that one can easily obtain the asymptotic stationary distribution P σ z (x, t → ∞).Indeed, from a careful inspection of Eq. ( 22) one can argue that for large time the contribution coming from I(s) is bounded by e −γt t 0 ds|I(s)|, thus decaying exponentially for large time.The only term that survives and contributes to the asymptotic stationary distribution is the one depending on e γs .Therefore, for any finite γ, the stationary distribution can be exactly evaluated as ds δ(x − σ cos(Ωs))e −γs = γ e 2 γ arcsin(x)/Ω + 1 e γ arccos(x)/Ω 2Ω (e πγ/Ω − 1) which, as expected, is an even function of x.Interestingly, all (non-vanishing) stationary moments can be easily computed from the integral representation of the probability, namely while the odd moments are identically vanishing.In particular, when γ → ∞, x 2n P σ z → 1 for all n, confirming the fact that the stationary distribution converges to [δ(x − 1) + δ(x + 1)]/2.

Hopping particle
A free hopping quantum particle propagates in a lattice with a ballistic spreading.However, there are ways to prevent or slow down the propagation as, for instance, adding a disorder potential which induces Anderson localization [60,61].Here, we show that the quantum Zeno effect due to the coupling of the hopping particle to a measurement apparatus can also results into a slowdown of the particle propagation [62,63].Related protocols have been studied in the context of quantum stochastic resetting, in which the hopping particle is reset to the initial state with a certain probability [64].
We consider a simple hopping fermion on a 1D lattice, whose hamiltonian reads with periodic boundary conditions (i.e.c j+L = c j ).Here the lattice dimension L (even) plays the role of an infrared cutoff.In the following we will take the limit L → ∞ whenever it will be unambiguous.The Fermionic operator satisfy the canonical anticommutation relations {c i , c † j } = δ ij .We define the Fourier modes operators such that {c p , c † q } = δ pq , and k ∈ {−π, −π + 2π/L, −π + 4π/L, . . ., π}.The Hamiltonian become diagonal in the Fourier representation, i.e.H = k ε(k)c † k ck , where ε(k) = −2Ω cos k.We now restrict the problem to the single-particle sector of the hamiltonian, we can define the states |j = c † j |∅ and | k = c † k |∅ , which represents the particle in position j or with momentum k respectively.Notice that j| k = exp(−ikj)/ √ L is the normalised wave function.Since H commute with he total number of particles, the unitary dynamics can be restricted in such sector and it is governed by the following single-particle Hamiltonian We consider the particle initial localised at the origin j = 0 of our lattice, namely |ψ ξ (0) = |0 for all trajectories ξ.We then suppose to continuously measure, with a rate γ, the position operator q = j j |j j| .
We are thus interested in the displacement of the particle along each single trajectory, however, for symmetry reasons, when no measurement occurs (i.e.γ = 0) the probability function of the outcome of q t is time independent, namely δ (x − q t ) = δ(0).This is not in contradiction with the expected ballistic spreading under the free evolution, which can be extracted when observing even power of q (see Appendix C).Notice that, this is not true anymore when γ = 0. We are thus interested in the probability distribution function of the particle displacement itself, namely where in this case the no-click contribution is trivially given by P (0) q (x; t) = e −γt δ(x), and we used the following results for the evolution amplitudes in the thermodynamic limit with J n (x) being the Bessel function of the first kind.Now the transition probability matrix reads T i,j (t) = J 2 i−j (2Ωt), which is a circulant matrix due to translational invariance.Let us define the not normalised eigenvectors whose component are where we identified the eigenvalues of the transfer matrix as with ω k = 4Ω sin(k/2).Following Section 2, we can easily solve the integral equation (13) for and we can identify with L [I n,k ](z) the terms of the series expansion, thus finally getting With a simple generalization of equation (10), we can write the click contribution to the probability distribution as Upper panels: the time-dependent probability distribution P q (x; t) for an hopping particle (Ω = 1).Red lines represents the standard deviation ± x 2 1/2 Pq as in Eq. (42).We set γ = 0.25 (left), γ = 1.0 (center) and γ = 5.0 (right).Lower panels: the expectation value ψ ξ (t)|q|ψ ξ (t) for 20 trajectories and the same values of γ.
where, thanks to the properties of the Bessel functions, the delta contribution reduces to D(x, j, t − s) = δ(x − j), which basically says that the probability has only support on x ∈ Z.The entire probability distribution P q (x; t) is represented in Figure 4 together with some representative quantum trajectories for different values of γ.
We have previously observed that the first moment of the probability distribution is identically vanishing due to the inversion symmetry.The second moment instead gets a nontrivial contribution from the P (n>0) q part, thus reading which can be further simplified using the identity j j 2 e −ijk = − j ∂ 2 k e ijk = −2πδ (k), leading to The second moment does behave differently depending on the time-scale, showing the following scaling behaviour Let us mention that the fluctuations of the P q (x; t) distribution cannot give information about the ballistic behavior at γ = 0, due to the inversion symmetry.As a matter of fact, what Eq. (43) gives us is the leading term for γ = 0.In the asymptotic regime γt 1, it does coincide (as expected) with x P q 2 confirming the diffusive behavior for any finite measurement rate.However, in the γt 1 regime, the O(γ 0 ) term is missing, and it is recovered in the expansion of x P q 2 (see Appendix C).

Conclusions
In this article, we have presented an approach for studying the dynamics of quantum systems undergoing unitary evolution and continuous monitoring.Our approach goes beyond the traditional Lindblad master equation and provides a more complete picture of the system's behavior by considering the entire ensemble of stochastic quantum trajectories.
Specifically, we have developed an analytical tool to compute the probability distribution of the expectation value of a given observable over the ensemble of quantum trajectories.We obtained exact formulas to evaluate this probability distribution and its moments for two paradigmatic examples: a single qubit subjected to magnetization measurements, and a free hopping particle subjected to position measurements.
Our results demonstrate that the probability distribution of expectation values can exhibit non-trivial features that are missed by traditional approaches, highlighting the full properties contained in the set of quantum trajectories in the analysis of continuously monitored quantum systems.
Note added.-Whilecompleting this work, the preprint [65] appeared, dealing with topics related to the ones discussed by us.
Combining the previous expansion with the unitary part of the evolution, and taking the continuum limit dt → 0 with γ fixed, we finally get the following Lindblad master equation with H = −Jσ x .This equation can be easily solved by expanding the density operator in the basis of Pauli matrices where m α (t) = Tr[σ α ρ(t)] and I is the 2 × 2 identity matrix.The Lindblad equation becomes a linear differential equation for the three components of the magnetisation (m x , m y , m z ) T , which reads where Ω = 2J.The z component evolves following the differential equation of a damped harmonic oscillator In the case of the hopping particle, to analyze the dynamical map averaged over quantum trajectories, we can again reframe the measurement procedure as a Lindblad equation for the averaged density matrix ρ(t).When we perform projective measurements of the particle's position at a rate γ, the average state ρ undergoes the following transformation in accordance with the usual rules of quantum mechanics ρ(t) → (1 − γdt) ρ(t) + γdt j π j ρ(t)π j . (C.1) Here, π j = |j j| represents the projectors over the lattice sites.The Lindblad master equation can be obtained by taking the limit dt → 0 with a fixed value of γ, resulting in Now, we collect few results for the first moment of some relevant observables.In particular, from the probability distribution of a generic power q m , by using the generic equation (18) we explicitly get x P q m = Tr{ρ(t)q m } = e −γt π −π dk 2π j j m e −ijk I k (t) .
(C.4) By exploiting j j m e −ijk = i m 2πδ (m) (k), we obtain As expected, x Pq = 0, while the first non-vanishing average occurs for m = 2, giving x P q 2 = 4Ω 2 (γt + e −γt − 1) thus, we recover the ballistic behavior at early time (or for γ = 0), whilst a diffusive behavior for any γ = 0 at large time, with a diffusion constant D γ = 2Ω 2 /γ such that x q 2 ∼ 2D γ t.This behavior has been observed in the many-body description of the model, in particular studying the particle current after quenching an initial domain wall configuration [24].Due to its linear nature, this exact result can be compared with the solution of the Lindblad equation.Indeed, the average of q 2 is also given by x P q 2 = j j 2 n(j, t), where n(j, t) = Tr{ρ(t)π j } can be evaluated numerically (see Fig. C1 upper panels), or analytically taking the average over the probability distribution associated to π j P π j (x; t) = δ ( ψ ξ (t)|π j |ψ ξ (t) − x).notice that, as expected, j n(j, t) = 1, since I 0 (t) = e γt .