Continuously monitored quantum systems beyond Lindblad dynamics

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–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–12] and have led to the development of novel applications in fields such as quantum information processing [13, 14] and quantum thermodynamics [15–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 [24–50]. 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 [51, 52]. 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 [53, 54].

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 We consider the case where projective measurements occur at a fixed rate γ. This choice seems quite natural, and corresponds to a Poissonian statistics of waiting times between successive measurements. However, the technique we introduce could, in principle, be applied to other waiting time statistics. 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, 55–57]. 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.

Let us consider an N-level quantum system described by a time-independent Hamiltonian H, whose unitary evolution is governed by $U(t) = e^{-i H t}$. 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

$$\begin{align} A = \sum_{a = 1}^N \nu_{a}|{a}\rangle\langle{a}|, \end{align} \tag{ 1 }$$

such that $[A,H]\neq 0$. At time t = 0 the system has been prepared in a fixed state $|{\psi(0)}\rangle$ corresponding to an eigenstate $|{a_0}\rangle$ of A. We are interested in evaluating the probability distribution of the expectation value of a generic observable $O = \sum_a o_a |{a}\rangle\langle{a}|$ commuting with A, i.e.

$$\begin{align} P_O\left(x;t\right) = \overline{\delta\left(\left\langle \psi_\xi\left(t\right) | O | \psi_\xi\left(t\right) \right\rangle-x\right)} \, , \end{align} \tag{ 2 }$$

with t > 0. Here, the over-line is indicating the average over the quantum trajectories, where each trajectory comprises a sequence of unitary time evolutions interspersed with instantaneous projective measurements of A at random times, and is labeled by an integer ξ. A sketch of our dynamical protocol is shown in figure 1. Considering this protocol, 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,\dots,n$ such that $0 \lt s_1\lt\cdots\lt s_n\lt t$. We define

$$\begin{align} P_O\left(x;t\right) = P^{\left(n>0\right)}_O\left(x;t\right) + P^{\left(0\right)}_O\left(x;t\right) \, , \end{align} \tag{ 3 }$$

where $P^{(0)}$, $P^{(n\gt0)}$ are respectively the no-click (n = 0) and click (n > 0) contributions to the probability distribution of $\left\langle \psi_\xi(t) | O | \psi_\xi(t) \right\rangle$. Now we proceed with the rewriting of these two quantities, explicitly specifying for $P^{(n\gt0)}$ the contribution coming from events with n measurements and summing over $n = 1,2,3,\dots$ (see [58] for a similar approach). Moreover, we have to integrate over the possible measurement times s_j and summing over their possible outcomes a_j . We obtain

$$\begin{align} \begin{cases} P^{\left(0\right)}_O\left(x;t\right) = e^{-\gamma t}\mathscr{D} \left(x,a_0,t\right) \\ \quad \\ P^{\left(n>0\right)}_O\left(x;t\right) = \sum_{n = 1}^\infty \gamma^n e^{-\gamma t} \sum_{a_1,\dots, a_n} \mathcal{T} \int_0^{t}\mathrm{d} s_1\dots \mathrm{d} s_n \, \, p\left(x,t| s_n,a_n; \dots; s_1,a_1;0,a_0\right), \end{cases} \end{align} \tag{ 4 }$$

where

$$\begin{align} \mathscr{D} \left(x,a,t\right)\equiv\delta \left(x-\left\langle a | {U}^{\dagger}ger\left(t\right) {O}{U}\left(t\right) | a \right\rangle \right) \end{align} \tag{ 5 }$$

and $p(x,t| s_n,a_n; {\ldots}; 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 \rangle \}$ of A at times $\{s_j\}$. $\mathcal{T}$ denotes the time-ordered product. The term $\gamma^n e^{-\gamma 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 labeled by $1,2{\ldots} \, n$ occurs at fixed times $s_1 \lt s_2 \lt {\cdots} \, \lt s_n$). Since the measurement process is Markovian, we have that

Zoom In Zoom Out Reset image size

$$\begin{align} p\left(x,t| s_n,a_n; \dots; s_1,a_1;0,a_0\right) = \mathscr{D}\left(x,a_n,t-s_n\right) p\left(s_n,a_n|s_{n-1},a_{n-1}\right) \dots p\left(s_1,a_1|0,a_0\right), \end{align} \tag{ 6 }$$

where the delta contribution can be rewritten as

$$\begin{align} \mathscr{D} \left(x,a,t\right) = \delta \left(x - \sum_{a^{^{\prime}} = 1}^{N} o_{a^{^{\prime}}} \, p\left(t,a^{^{\prime}}|0,a\right) \right), \end{align} \tag{ 7 }$$

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 ³

$$\begin{align} p\left(s^{^{\prime}},a^{^{\prime}}|s,a\right) \equiv \mathbb{T}_{a^{^{\prime}},a}\left(s^{^{\prime}}-s\right) = |{\left\langle a^{^{\prime}} | U\left(s^{^{\prime}}-s\right) | a \right\rangle}|^2. \end{align} \tag{ 8 }$$

Let us assume that the transition matrix $\mathbb{T}$ can be diagonalized, i.e. $\mathbb{T}(s^{^{\prime}}-s) = \mathbb{V} \, \mathbb{D}(s^{^{\prime}}-s) \,\mathbb{V}^{\dagger}$ with eigenvalues $\mathbb{D}_{\alpha,\beta}(t) = d_\alpha(t)\delta_{\alpha,\beta}$, and time-independent unitary matrix $\mathbb{V}$ whose columns are the orthonormal eigenvectors ⁴ . Notice that, since $\mathbb{T}(t)$ is a bistochastic matrix, $d_1(t) = 1$ is the largest eigenvalue for any time t, meanwhile as a consequence of the Perron–Frobenius theorem $|{d_\alpha (t)}| \unicode{x2A7D} 1$ for $\alpha = 2,\dots,N$, where the equality holds in case of degeneracy [60, 61]. In addition, since $\sum_a \mathbb{T}_{a,b} = \sum_b \mathbb{T}_{a,b} = 1$, the eigenvector associated to largest eigenvalue is simply given by $\mathbb{V}_{a,1} = 1/\sqrt{N}$. Now let us focus on the nontrivial $P^{(n\gt0)}_O(x;t)$. For each order n and fixed states $\{a_0, a_n\}$, the sum over all the possible intermediate outcomes can be rewritten as

$$\begin{align} \sum_{a_1,\dots,a_{n-1}} \mathbb{T}_{a_{n},a_{n-1}}\left(s_n-s_{n-1}\right) \cdots \mathbb{T}_{a_1,a_0}\left(s_1\right) = \sum_{\alpha = 1}^{N}\mathbb{V}_{a_n,\alpha} \prod_{j = 1}^{n} d_{\alpha}\left(s_{j}-s_{j-1}\right) \left(\mathbb{V}^{\dagger}\right)_{\alpha,a_0}, \end{align} \tag{ 9 }$$

with starting time $s_0 = 0$. We thus get

$$\begin{align} P^{\left(n>0\right)}_O\left(x;t\right) = e^{-\gamma t}\sum_{n = 1}^\infty \left[ \int_0^t \gamma \mathrm{d} s_n \sum_{a_n = 1}^{N} \mathscr{D} \left(x,a_n,t-s_n\right) \sum_{\alpha = 1}^N \mathbb{V}_{a_n,\alpha} I_{n,\alpha}\left(s_n\right) \left(\mathbb{V}^{\dagger}\right)_{\alpha,a_0} \right] \end{align} \tag{ 10 }$$

where we defined the time-ordered integrals

$$\begin{align} I_{n,\alpha}\left(s_n\right) = \int_0^{s_n} \gamma \mathrm{d} s_{n-1} \, d_\alpha\left(s_n-s_{n-1}\right) \, \dots \, \int_0^{s_2}\gamma \mathrm{d} s_1 \, d_\alpha\left(s_2-s_1\right)d_\alpha\left(s_1\right), \end{align} \tag{ 11 }$$

which satisfy the following recursive relations

$$\begin{align}\left\{ \begin{array}{@{}l} I_{1,\alpha}\left(s_1\right) = d_\alpha\left(s_1\right),\\[3pt] I_{n,\alpha}\left(s_n\right) = \displaystyle\int_0^{s_n}\gamma \mathrm{d} s_{n-1} \, d_\alpha\left(s_n-s_{n-1}\right) I_{n-1,\alpha}\left(s_{n-1}\right), \quad n>1 . \end{array}\right. \end{align} \tag{ 12 }$$

In equation (10) the sum over n can be safely taken, since s_n and a_n are dummies variables. Thus, we can introduce $I_{\alpha}(s) = \sum_{n = 1}^{\infty} I_{n,\alpha}(s)$ which satisfy the following linear Volterra integral equation of the second kind

$$\begin{align} I_\alpha\left(s\right) - d_\alpha\left(s\right) = \int_0^s \gamma \mathrm{d} s^{^{\prime}} \, d_\alpha\left(s-s^{^{\prime}}\right)I_\alpha\left(s^{^{\prime}}\right). \end{align} \tag{ 13 }$$

Exploiting the Laplace transform of the integral kernel, i.e. $\mathscr{L}[d_\alpha](z) = \int_0^\infty e^{-zt} d_\alpha(t) \mathrm{d} t$, we can easily solve the previous equation obtaining

$$\begin{align} I_\alpha\left(s\right) = \mathscr{L}^{-1}\left[ \frac{\mathscr{L}\left[d_\alpha\right]\left(z\right)}{1-\gamma\mathscr{L}\left[d_\alpha\right]\left(z\right)}\right]\left(s\right), \end{align} \tag{ 14 }$$

where $\mathscr{L}^{-1}[\dots]$ is the inverse Laplace transform. As expected, it easy to show that

$$\begin{align} \mathscr{L}\left[I_\alpha\right]\left(z\right) = \frac{\mathscr{L}\left[d_\alpha\right]\left(z\right)}{1-\gamma\mathscr{L}\left[d_\alpha\right]\left(z\right)} = \sum_{n = 1}^{\infty} \gamma^{n-1}\mathscr{L}\left[d_\alpha\right]^n\left(z\right) = \sum_{n = 1}^{\infty}\mathscr{L}\left[I_{n,\alpha}\right]\left(z\right), \end{align} \tag{ 15 }$$

thus identifying $\mathscr{L}[I_{n,\alpha}](z) = \gamma^{n-1}\mathscr{L}[d_\alpha]^n(z)$. We can finally collect all those results, further simplify and finding the following general expression for the probability distribution

$$\begin{align} P^{\left(n>0\right)}_O\left(x;t\right) = e^{-\gamma t} \sum_{a = 1}^{N} \sum_{\alpha = 1}^N \int_0^t \gamma \mathrm{d} s \, \mathscr{D} \left(x,a,t-s\right) \mathbb{V}_{a,\alpha} I_{\alpha}\left(s\right)\left(\mathbb{V}^{\dagger}\right)_{\alpha,a_0}. \end{align} \tag{ 16 }$$

Let us finally mention that using this expression we can easily compute the kth moment of the distribution, defined as $\left\langle x^{k} \right\rangle_{P_{O}} = \int \mathrm{d} x \, x^k P_O (x;t) = \overline{\left\langle \psi_\xi(t) | O | \psi_\xi(t) \right\rangle^k}$. The first moment (k = 1) is particularly simple, since it does reduce to the expectation value of O over the averaged density matrix $\overline{| \psi_{\xi}(t) \rangle \langle{\psi_{\xi}(t)}|}$ whose dynamics is fully described by the Lindblad equation (see appendix C). In particular, using that $\sum_{b}\mathbb{T}_{a,b}(t-s) \mathbb{V}_{b,\alpha} =$ $d_{\alpha}(t-s) \mathbb{V}_{a,\alpha}$, we easily get

$$\begin{align} \left\langle x^{} \right\rangle_{P_{O}} = e^{-\gamma t} \sum_{a = 1}^{N} o_{a} \, \mathbb{T}_{a,a_0}\left(t\right) + e^{-\gamma t} \sum_{a = 1}^{N} o_a \sum_{\alpha = 1}^{N} \mathbb{V}_{a,\alpha} \left[\int_0^t \gamma \mathrm{d} s \, d_\alpha\left(t-s\right) I_{\alpha}\left(s\right) \right] \left(\mathbb{V}^{\dagger}\right)_{\alpha,a_0}, \end{align} \tag{ 17 }$$

which can be further simplified exploiting equation (13), finally obtaining

$$\begin{align} \left\langle x^{} \right\rangle_{P_{O}} = e^{-\gamma t} \sum_{a = 1}^{N} o_a \sum_{\alpha = 1}^{N} \mathbb{V}_{a,\alpha} I_{\alpha}\left(t\right) \left(\mathbb{V}^{\dagger}\right)_{\alpha,a_0}. \end{align} \tag{ 18 }$$

Let us stress that the simplified expression in equation (18) applies only for averages of operators O which are diagonal in the eigenbasis of the monitored observable A (i.e. when $[O,A] = 0$). 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 equation (16). Finally, due to the properties of the transfer matrix $\mathbb{T}$, the stationary limit $t\to\infty$ of the previous average can be easily taken; in fact, only α = 1, with $I_{1}(t) = e^{\gamma t}$, will contribute to the sum over α, leading to $\left\langle x^{} \right\rangle_{P_{O}} \to \frac{1}{N} \sum_{a = 1}^{N} o_a$, where we used the fact that $\mathbb{V}_{a,1} = 1/\sqrt{N}$. As expected from the Lindblad dynamics, it does correspond to the expectation value of the operator O over the infinite temperature state $\frac{1}{N} \sum_{a = 1}^{N} | a \rangle\langle{a}|$.

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

$$\begin{align} {H} = -J \sigma^x , \end{align} \tag{ 19 }$$

where σ^α are the Pauli matrices with $\alpha = x,y,z$ and $[{\sigma^\alpha}, {\sigma^\beta}] = 2i\epsilon_{\alpha\beta\gamma}\sigma^\gamma$. The system is continuously monitored along z with a rate γ (therefore ${A} = \sigma^z$). We will denote with $| \sigma \rangle = | +1 \rangle,| -1 \rangle$ the eigenstates of σ^z with eigenvalue $\sigma = \pm 1$ respectively. We consider a two-level spin which starts from $| +1 \rangle$ and we want to evaluate

$$\begin{align} P_{\sigma^z}\left(x;t\right) = \overline{\delta\left(\left\langle \psi_\xi\left(t\right) | \sigma^z | \psi_\xi\left(t\right) \right\rangle-x\right)} = P^{\left(n>0\right)}_{\sigma^z}\left(x;t\right)+P^{\left(0\right)}_{\sigma^z}\left(x;t\right) \, , \end{align} \tag{ 20 }$$

i.e. the probability distribution of the expectation value of ${O} = \sigma^z$. It is easy to find that $\mathbb{T}(t) = \text{cos}(J t)^2 \mathbb{I} + \text{sin}(J t)^2 \sigma^x$, and $\mathbb{V} = (\sigma^z +\sigma^x)/\sqrt{2}$ whose columns are the eigenvectors of σ^x . As expected $d_1 = 1$ and $d_2(t) = \text{cos}(\Omega t)$, where $\Omega = 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 equation (12), then, since $d_1 = 1$ we immediately get $I_1(s) = e^{\gamma s}$. In the following, we will simplify the notation $I_2(s) = I(s)$. Finally, we notice that $\mathscr{D} (x,\sigma_n,t-s_n) = \delta \big(x-\sigma_n\text{cos}(\Omega(t-s_n)) \big)$, from which we can compute the no-click contribution

$$\begin{align} P^{\left(0\right)}_{\sigma^z}\left(x;t\right) = e^{-\gamma t} \, \delta \left(x-\text{cos}\left(\Omega t\right) \right). \end{align} \tag{ 21 }$$

In order to find the full probability distribution, we can use equation (16) to write

$$\begin{align} P^{\left(n>0\right)}_{\sigma^z}\left(x;t\right) = \frac{e^{-\gamma t}}{2} \int_0^t \gamma \mathrm{d} s \sum_{\sigma = \pm 1} \delta\left(x-{\sigma}\text{cos}\left(\Omega\left(t-s\right)\right)\right) \left[ e^{\gamma s} + \sigma I\left(s\right) \right]. \end{align} \tag{ 22 }$$

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)

$$\begin{align} \ddot{I} = \gamma \dot{I} - \Omega^2 I \end{align} \tag{ 23 }$$

with the initial conditions $I(0) = 1$ and $\dot{I}(0) = \gamma$. This is the equation of motion of an 'anti-damped' harmonic oscillator (since γ > 0). The solution can be found easily using standard methods [62], obtaining

$$\begin{align} I\left(s\right) = e^{\gamma s/2}\left[\text{cosh}\left(\frac{\Gamma s}{2}\right)+\frac{\gamma}{\Gamma}\text{sinh}\left(\frac{\Gamma s}{2}\right)\right], \end{align} \tag{ 24 }$$

where we defined the following characteristic frequency $\Gamma = \sqrt{\gamma ^2-4 \Omega^2}$. As expected, we find three different regimes depending on the sign of $\gamma^2-4\Omega^2$. Let us analyze the behavior of the function $e^{-\gamma s/2}I(s)$. When $\gamma\lt2\Omega$, we have that Γ is an imaginary quantity. Therefore, we find a regime in which the function oscillates with a frequency of $\Omega\sqrt{1-\gamma^2/(4\Omega^2)}$, which is lower than the natural frequency Ω. When $\gamma = 2\Omega$, we get the critical regime for which the function $e^{-\gamma s/2}I(s)$ becomes linear, $\left(1+\gamma s/2\right)$. Finally, if $\gamma\gt2\Omega$, it grows exponentially with a rate $\gamma/2\sqrt{1-4\Omega^2/\gamma^2}$ smaller than $\gamma/2$.

Let us finally consider the Dirac-delta function contribution. Solving for $x-\sigma \text{cos}(\Omega(t-s)) = 0$, we easily get

$$\begin{align} \tilde{s}_k = t \pm \frac{1}{\Omega}\text{arccos}\left(x\sigma\right)+\frac{2\pi k}{\Omega},\quad k \in \mathbb{Z}. \end{align} \tag{ 25 }$$

Using the properties of the Dirac-delta distribution, we can rewrite it as

$$\begin{align} \delta\left(x-\sigma \text{cos}\left(\Omega\left(t-s\right)\right)\right) = \sum_{k \in \mathbb{Z}} \frac{\delta\left(s-\tilde{s}_k\right)}{\Omega\sqrt{1-x^2}}. \end{align} \tag{ 26 }$$

Since in equation (22), the integral is taken over a finite interval $[0,t]$, we need to constraint the solutions $\tilde s_k$ in equation (25) only to those falling in that region. However, we may relax that condition by explicitly introducing the Heaviside step function $\theta(x)$, such that equation (22) finally reads

$$\begin{align} P^{\left(n>0\right)}_{\sigma^z}\left(x;t\right) = \frac{\gamma e^{-\gamma t}}{2\Omega \sqrt{1-x^2}} \sum_{\sigma = \pm 1} \sum_{k \in \mathbb{Z}} \left[ e^{\gamma \tilde{s}_k} + \sigma I\left(\tilde{s}_k\right) \right]\theta\left(t-\tilde{s}_k\right) \theta\left(\tilde{s}_k\right) \, , \end{align} \tag{ 27 }$$

where $\tilde s_k$ implicitly depends on x and t as well. In figure 2, we plot the entire distribution $P_{\sigma^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 localized at $x = \text{cos}(\Omega t)$ but its weight is exponentially suppressed in time. The click contribution is also showing a nontrivial asymmetric evolution and has a discontinuity corresponding to the value of the no-click delta peak (see lower panels). Its behavior strongly depends on γ indeed: in the oscillatory regime ($\gamma \lt 2\Omega$), the probability distribution function bounces back and forth between the two extremes $x = \pm 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 $\tau = 2/\gamma$. For $\gamma \gt 2\Omega$ the probability is not oscillating anymore and the relaxation time becomes $\tau = 2/(\gamma - \Gamma)$. Interestingly, as the measurement rate is getting higher, the time needs for the magnetization statistics to reach the equilibrium becomes larger and larger, diverging as $\tau \sim \gamma/\Omega^2$.

Zoom In Zoom Out Reset image size

The first moment of $P_{\sigma^z}(x;t)$, namely the magnetization average $\left\langle x^{} \right\rangle_{P_{{\sigma^z}}} = e^{-\gamma 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 equations (21) and (22), we can easily get the second moment of the distribution

$$\begin{align} \left\langle x^{2} \right\rangle_{P_{{\sigma^z}}} = e^{-\gamma t}\int_0^t \gamma \mathrm{d} s \text{cos}\left(\Omega\left(t-s\right)\right)^2 e^{\gamma s} + e^{-\gamma t}\text{cos}\left(\Omega t\right)^2, \end{align} \tag{ 28 }$$

which simplifies to

$$\begin{align} \left\langle x^{2} \right\rangle_{P_{{\sigma^z}}} = \frac{\gamma^2+2\Omega^2}{\gamma^2+4\Omega^2}+\frac{e^{-\gamma t}}{2\left(\gamma^2+4\Omega^2\right)}\left[2\gamma\Omega\text{sin}\left(2\Omega t\right)+\Omega^2\text{cos}\left(2\Omega t\right)\right], \end{align} \tag{ 29 }$$

and which cannot be computed using the a Lindblad approach since it corresponds to $\overline{\langle\psi_{\xi}(t)|\sigma_z |\psi_{\xi}(t) \rangle^2}$. In figure 3, we represents $\left\langle x^{} \right\rangle_{P_{{\sigma^z}}}, \left\langle x^{2} \right\rangle_{P_{{\sigma^z}}}$ as a function of time for different values of the measurement rate γ.

Zoom In Zoom Out Reset image size

Finally, let us mention that one can easily obtain the asymptotic stationary distribution $P_{\sigma^z}(x,t\to\infty)$. Indeed, from a careful inspection of equation (22) one can argue that for large time the contribution coming from I(s) is bounded by $e^{-\gamma t}\int_{0}^{t} \mathrm{d} s | 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^{\gamma s}$. Therefore, for any finite γ, the stationary distribution can be exactly evaluated as

$$\begin{align*} P_{\sigma^z}\left(x\right) = \frac{\gamma}{2}\sum_{\sigma = \pm 1}\int_0^\infty\mathrm{d} s\, \delta\left(x-\sigma\text{cos}\left(\Omega s\right)\right)e^{-\gamma s} = \frac{\gamma \left[e^{2\,\gamma \text{arcsin}\left(x\right)/\Omega}+1\right] e^{\gamma \text{arccos}\left(x\right)/\Omega}}{2 \Omega \left(e^{\pi \gamma/\Omega}-1\right) \sqrt{1-x^2}}, \end{align*}$$

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

$$\begin{align} \left\langle x^{2n} \right\rangle_{P_{{\sigma^z}}} = \int_{0}^{\infty} ds \, e^{-s} \text{cos}^{2n}\left(s \Omega/\gamma\right) = \frac{_2 F_1\left(-2n,-n-i \gamma/\left(2\Omega\right);1-n-i \gamma/\left(2\Omega\right);-1\right)}{2^{2n} \left(1 - 2 i n \Omega/\gamma\right)}, \end{align} \tag{ 30 }$$

while the odd moments are identically vanishing. In particular, when $\gamma\to\infty$, $\left\langle x^{2n} \right\rangle_{P_{{\sigma^z}}}\to 1$ for all n, confirming the fact that the stationary distribution converges to $[\delta(x-1)+\delta(x+1)]/2$.

Observe that, having already taken the limit $t \rightarrow \infty$ to derive the stationary distribution in equation (3), the limit $\gamma \rightarrow \infty$ does not result in the distribution linked to Zeno's effect, namely $\delta(x-1)$. To clarify, the order of limits matters, and to obtain the Zeno's effect distribution, one have to first take the limit $\gamma \rightarrow \infty$.

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 [63, 64]. 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 [65, 66]. 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 [67].

We consider a simple hopping fermion on a 1D lattice, whose Hamiltonian reads

$$\begin{align} H = -\Omega \sum_{j = 1}^{L} \left( c^{\dagger}_j c_{j+1} + c^{\dagger}_{j+1} c_j \right) \end{align} \tag{ 31 }$$

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\to\infty$ whenever it will be unambiguous. The Fermionic operator satisfy the canonical anti-commutation relations $\{c_i,c^{\dagger}_j\} = \delta_{ij}$. We define the Fourier modes operators

$$\begin{align} \tilde c_{k} = \frac{1}{\sqrt{L}}\sum_j e^{ik j} c_j, \quad c_j = \frac{1}{\sqrt{L}} \sum_k e^{-ikj} \tilde c_k, \end{align} \tag{ 32 }$$

such that $\{\tilde c_p,\tilde c^{\dagger}_q\} = \delta_{pq}$, and $k \in \{-\pi,-\pi+2\pi/L,-\pi+4\pi/L,\dots, \pi\}$. The Hamiltonian become diagonal in the Fourier representation, i.e. $H = \sum_k \varepsilon(k) \tilde c^{\dagger}_k \tilde c_k$, where $\varepsilon(k) = -2 \Omega \text{cos} k$. We now restrict the problem to the single-particle sector of the Hamiltonian, we can define the states $| j \rangle = c^{\dagger}_j| \emptyset \rangle$ and $| \tilde k \rangle = \tilde c^{\dagger}_k| \emptyset \rangle$, which represents the particle in position j or with momentum k respectively. Notice that $\langle{j}|{\tilde k}\rangle = \text{exp}(-ikj)/\sqrt{L}$ is the normalized wave function. Since H commute with the total number of particles, the unitary dynamics can be restricted in such sector and it is governed by the following single-particle Hamiltonian

$$\begin{align} H = -\Omega\sum_{j = 1}^{L} \left( |{j}\rangle\langle{j+1}| + |{j+1}\rangle\langle{j}| \right) = \sum_k \varepsilon\left(k\right)| \tilde k \rangle\langle \tilde k |.\end{align} \tag{ 33 }$$

We consider the particle initial localized at the origin j = 0 of our lattice, namely $| \psi_{\xi}(0) \rangle = | 0 \rangle$ for all trajectories ξ. We then supposed to continuously measure, with a rate γ, the position operator

$$\begin{align} q = \sum_{j} j | j \rangle\langle j |. \end{align} \tag{ 34 }$$

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 $\langle q\rangle_t$ is time independent, namely $\overline{\delta\left(x-\langle q\rangle_t\right)} = \delta(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

$$\begin{align} P_q\left(x;t\right) = \overline{\delta\left(\left\langle \psi_\xi\left(t\right) | q | \psi_\xi\left(t\right) \right\rangle-x\right)} = P^{\left(n>0\right)}_q\left(x;t\right)+P^{\left(0\right)}_q\left(x;t\right), \end{align} \tag{ 35 }$$

where in this case the no-click contribution is trivially given by $P^{(0)}_q(x;t) = e^{-\gamma t} \delta(x)$, and we used the following results for the evolution amplitudes in the thermodynamic limit

$$\begin{align} \left\langle j | U\left(t\right) | l \right\rangle = \sum_{k} \left\langle j | e^{-i t H} | \tilde k \right\rangle\langle{\tilde k}|{l}\rangle & = \int_{-\pi}^\pi \frac{\mathrm{d} k}{2\pi} e^{i\left[k\left(j-l\right)-2t\Omega\text{cos} k\right]} = \left(-i\right)^{j-l}J_{j-l}\left(2\Omega t\right), \end{align} \tag{ 36 }$$

with $J_n(x)$ being the Bessel function of the first kind. Now the transition probability matrix reads $\mathbb{T}_{i,j}(t) = J^2_{i-j}(2\Omega t)$, which is a circulant matrix due to translational invariance. Let us define the not normalized eigenvectors whose component are $\mathbb{V}_{n,k} = e^{-ink}$, such that

$$\begin{align} \left(\mathbb{T}\left(t\right)\mathbb{V}\right)_{n,k} = \sum_j J^2_{n-j}\left(2\Omega t\right) e^{-ijk} = e^{-ink}\sum_l J^2_j\left(2\Omega t\right)e^{ilk} \equiv \mathbb{V}_{n,k} \, d_k\left(t\right) \end{align} \tag{ 37 }$$

where we identified the eigenvalues of the transfer matrix as $d_k(t) = \sum_l J^2_l(2\Omega t)e^{ikl} = J_{0}[\omega_k t]$, with $\omega_k = 4 \Omega \text{sin}(k/2)$. Following section 2, we can easily solve the integral equation (13) for $I_k(s) = \sum_{n = 1}^\infty I_{n,k}(s)$ with kernel $J_{0}[\omega_k t]$; the Laplace transform reads

$$\begin{align} \mathscr{L}\left[I_k\right]\left(z\right) = \frac{1}{\sqrt{z^2+\omega_k^2}-\gamma} = \sum_{n = 1}^{\infty}\frac{\gamma^{n-1}}{\left(z^2+\omega_k^2\right)^{n/2}}, \end{align} \tag{ 38 }$$

and we can identify with $\mathscr{L}[I_{n,k}](z)$ the terms of the series expansion, thus finally getting

$$\begin{align} I_{n,k}\left(s\right) = \sqrt{\pi } \gamma ^{n-1}\left(\frac{s}{2 |\omega_k|}\right)^{\frac{n-1}{2}}\frac{ J_{\frac{n-1}{2}}\left(s |\omega_k|\right)}{\Gamma \left(n/2\right)}. \end{align} \tag{ 39 }$$

With a simple generalization of equation (10), we can write the click contribution to the probability distribution as

$$\begin{align} P^{\left(n>0\right)}_{q}\left(x;t\right) = e^{-\gamma t} \sum_{j} \int_0^t \gamma \mathrm{d} s \int_{-\pi}^\pi \frac{\mathrm{d} k}{2\pi} \mathscr{D} \left(x,j,t-s\right) e^{-ijk} I_{k}\left(s\right), \end{align} \tag{ 40 }$$

where, thanks to the properties of the Bessel functions, the delta contribution reduces to $\mathscr{D} (x,j,t-s) = \delta(x-j)$, which basically says that the probability has only support on $x\in\mathbb{Z}$. The entire probability distribution $P_q(x;t)$ is represented in figure 4 together with some representative quantum trajectories for different values of γ.

Zoom In Zoom Out Reset image size

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\gt0)}_q$ part, thus reading

$$\begin{align} \left\langle x^{2} \right\rangle_{P_{q}} = e^{-\gamma t}\int_0^t \gamma \mathrm{d} s\, \int_{-\pi}^\pi \frac{\mathrm{d} k}{2\pi} \left[\sum_j j^2 e^{-ij k}\right]I_k\left(s\right), \end{align} \tag{ 41 }$$

which can be further simplified using the identity $\sum_j j^2 e^{-ijk} = -\sum_j \partial_k^2 e^{ijk} = - 2 \pi \delta^{^{\prime\prime}}(k),$ leading to

$$\begin{align} \left\langle x^{2} \right\rangle_{P_{q}} = -e^{-\gamma t }\int_0^t \gamma \mathrm{d} s \, \left[\partial^2_k I_k\left(s\right)\right]\bigr\rvert_{k = 0} = \frac{4\Omega^2}{\gamma^2} \left[\left(\gamma t-2\right)+e^{-\gamma t} \left(\gamma t+2\right)\right]. \end{align} \tag{ 42 }$$

The second moment does behave differently depending on the time-scale, showing the following scaling behavior

$$\begin{align} \left\langle x^{2} \right\rangle_{P_{q}} \sim \begin{cases} \gamma \Omega^2t^3\quad \gamma t \ll 1\\ \frac{4\Omega^2t}{\gamma} \quad \gamma t \gg 1 \end{cases}.\end{align} \tag{ 43 }$$

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 equation (43) gives us is the leading term for γ ≠ 0. In the asymptotic regime $\gamma t \gg 1$, it does coincide (as expected) with $\left\langle x^{} \right\rangle_{P_{q^2}}$ confirming the diffusive behavior for any finite measurement rate. However, in the $\gamma t \ll 1$ regime, the $O(\gamma^0)$ term is missing, and it is recovered in the expansion of $\left\langle x^{} \right\rangle_{P_{q^2}}$ (see appendix C).

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. However, it is worth noting that obtaining this probability distribution is challenging in any experimental setup. Indeed, like the entanglement entropy, it requires repetitive measurements along individual quantum trajectories, which brings up the well-known issue of post-selection. Future studies could explore extending our methodology to continuous weak measurement scenarios, which involve non-Hermitian dynamics, though it remains an open question how to effectively generalize our formalism to such cases.

This work was supported by the PNRR MUR Project PE0000023-NQSTI

Note added.—While completing this work, the preprint [68] appeared, dealing with topics related to the ones discussed by us.

All data that support the findings of this study are included within the article (and any supplementary files).

Here we collect some properties of the Bessel functions

$$\begin{align} J_n\left(t\right) = \frac{1}{2\pi} \int_{-\pi}^\pi \mathrm{d} z\, e^{i\left(nz -t\text{sin} z\right)}, \end{align} \tag{ A.1 }$$

which have been used all along the main text. Let's start by noticing that $J_{n}(t)^2$ are non-negative and normalized over $\mathbb{Z}$, i.e. $\sum_{n\in\mathbb{Z}}J_{n}(t)^{2} = 1$. In other words, they can be interpreted as probabilities over the infinite set of discrete events $n\in\mathbb{Z}$, with a parameter dependence t. They indeed satisfy the following property

$$\begin{align} \sum_{n\in\mathbb{Z}} x^{2n} J_{n}\left(t\right)^2 = J_{0}\left[i t \left(x-1/x\right)\right], \end{align} \tag{ A.2 }$$

from which we can easily define the generating function $F_t(\lambda)$ of the moments of the distribution, as

$$\begin{align} F_{t}\left(\lambda\right) \equiv \sum_{n\in\mathbb{Z}} e^{i\lambda n}J_{n}^{2}\left(t\right) = J_{0}\left[2t \text{sin}\left(\lambda/2\right)\right], \end{align} \tag{ A.3 }$$

such that $\sum_{n\in\mathbb{Z}} n^{k} J_{n}(t)^2 = (-i)^{k}\partial^{k}_{\lambda}F_t(\lambda)$. From the previous relation it is straightforward to show that

$$\begin{align} \int_{-\pi}^{\pi}\frac{\mathrm{d} \lambda}{2\pi} e^{-i \lambda n} J_{0}\left[2t \text{sin}\left(\lambda/2\right)\right] = J^2_n\left(t\right).\end{align} \tag{ A.4 }$$

Finally, another useful property of the Bessel functions is that it is possible to compute the Laplace Transform, which reads

$$\begin{align} \mathscr{L}\left[J_n\right]\left(z\right) = \frac{\left(z+\sqrt{z^2+1}\right)^{-n}}{\sqrt{z^2+1}}.\end{align} \tag{ A.5 }$$

When we are interested in the dynamical map averaged over the quantum trajectories, the measurement protocol outlined in the main text can be reformulated in terms of a Lindblad equation for the averaged density matrix $\rho(t) = \overline{| \psi_{\xi}(t) \rangle\langle \psi_{\xi}(t) |}$. In fact, averaging over different trajectories does correspond to relax both the information on whether the spin has been measured, and the result of the measurement itself. See [69] for some results of projective measurement-based dissipation descriptions of Lindblad equations for quantum spin systems.

In particular for a single spin undergoing projective measurements of σ^z (see section 3), the average state ρ transforms accordingly to

$$\begin{align} \rho\left(t\right) \to \left( 1-\frac{\gamma {\mathrm {d}}t}{2} \right) \rho\left(t\right) + \frac{\gamma {\mathrm {d}}t}{2} \sigma^{z} \, \rho\left(t\right) \, \sigma^{z}, \end{align} \tag{ B.1 }$$

where $\gamma {\mathrm {d}}t$ is the probability that the system is measured, after a discretization of the continuum time evolution has been applied.

Combining the previous expansion with the unitary part of the evolution, and taking the continuum limit ${\mathrm {d}}t\to0$ with γ fixed, we finally get the following Lindblad master equation

$$\begin{align} \partial_t \rho = - i \left[ H, \rho\right] + \frac{\gamma}{2} \left( \sigma^{z} \, \rho \, \sigma^{z} - \rho \right), \end{align} \tag{ B.2 }$$

with $H = -J \sigma^{x}$. This equation can be easily solved by expanding the density operator in the basis of Pauli matrices

$$\begin{align} \rho\left(t\right) = \frac{1}{2}\mathbb{I} + \frac{1}{2}\sum_{\alpha = 1}^3 m_{\alpha}\left(t\right) \sigma^{\alpha}, \end{align} \tag{ B.3 }$$

where $m_{\alpha}(t) = \mathrm{Tr}[\sigma^{\alpha}\rho(t)]$ and $\mathbb{I}$ is the $2\times 2$ identity matrix. The Lindblad equation becomes a linear differential equation for the three components of the magnetization $(m_x,m_y,m_z)^T$, which reads

$$\begin{align} \partial_t \begin{pmatrix} m_{x}\\ m_{y}\\ m_z \end{pmatrix} = \begin{pmatrix} -\gamma & 0 & 0\\ 0 & -\gamma & \Omega\\ 0 & -\Omega & 0 \end{pmatrix} \begin{pmatrix} m_{x}\\ m_{y}\\ m_z \end{pmatrix} \end{align} \tag{ B.4 }$$

where $\Omega = 2J$. The z component evolves following the differential equation of a damped harmonic oscillator $\partial^2_t m_z = - \Omega^2 m_z -\gamma \partial_t m_z$, with initial condition $m_z(0) = 1$, $\partial_t m_z(0) = 0$. From the solution of such equation we easily recover the result of the main text, namely

$$\begin{align} m_z\left(t\right) = e^{-\gamma t /2} \left[\text{cosh}\left(\frac{\Gamma t}{2}\right)+\frac{\gamma}{\Gamma}\text{sinh}\left(\frac{\Gamma t}{2}\right)\right] \equiv \left\langle x^{} \right\rangle_{P_{\sigma^z}} \end{align} \tag{ B.5 }$$

with $\Gamma = \sqrt{\gamma ^2-4 \Omega^2 }$. In addition, we also gets

$$\begin{align} m_x\left(t\right) = 0, \quad m_y\left(t\right) = e^{-\gamma t /2} \frac{2\Omega}{\Gamma}\text{sinh}\left(\frac{\Gamma t}{2}\right). \end{align} \tag{ B.6 }$$

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 $\rho(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

$$\begin{align} \rho\left(t\right) \to \left(1-\gamma \mathrm{d} t\right)\rho\left(t\right)+\gamma \mathrm{d} t\sum_j \pi_j \rho\left(t\right) \pi_j \, . \end{align} \tag{ C.1 }$$

Here, $\pi_j = | j \rangle\langle j |$ represents the projectors over the lattice sites. The Lindblad master equation can be obtained by taking the limit $\mathrm{d} t \to 0$ with a fixed value of γ, resulting in

$$\begin{align} \partial_t \rho = -i \left[H,\rho\right]+\gamma\left[\sum_j \pi_j \rho \pi_j - \rho\right], \end{align} \tag{ C.2 }$$

with H as in equation (33). In the case of a system with finite size, we can solve numerically the dynamical map

$$\begin{align} \rho\left(t+\mathrm{d} t\right) = \rho\left(t\right)-i\mathrm{d} t \left(H\rho-\rho H\right)+\gamma \mathrm{d} t \left[\sum_j \pi_j \rho\left(t\right) \pi_j -\rho\left(t\right)\right], \end{align} \tag{ C.3 }$$

and define the site densities $n(j,t) = \mathrm{Tr}\{{\rho(t)\pi_j}\}$ as depicted in the upper four panels in figure C1. In accordance to what have been discussed in the main text, the averages over the probability distribution coincides with the expectation values taken over the averaged state.

Zoom In Zoom Out Reset image size

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

$$\begin{align} \left\langle x^{} \right\rangle_{P_{q^m}} = \mathrm{Tr}\{{\rho\left(t\right)q^m}\} = e^{-\gamma t} \int_{-\pi}^{\pi} \frac{\mathrm{d} k}{2\pi} \sum_{j} j^m e^{-i j k} I_{k}\left(t\right) \, . \end{align} \tag{ C.4 }$$

By exploiting $\sum_{j} j^m e^{-i j k} = i^m 2\pi \delta^{(m)}(k)$, we obtain

$$\begin{align} \left\langle x^{} \right\rangle_{P_{q^m}} = i^m e^{-\gamma t} \partial^{m}_{k} I_k\left(t\right)|_{k = 0}. \end{align} \tag{ C.5 }$$

As expected, $\left\langle x^{} \right\rangle_{P_{q}} = 0$, while the first non-vanishing average occurs for m = 2, giving

$$\begin{align} \left\langle x^{} \right\rangle_{P_{q^2}} = \frac{4 \Omega^2 \left(\gamma t+e^{-\gamma t}-1\right)}{\gamma ^2} \sim \begin{cases} 2 \Omega^2 t^2 \quad \gamma t \ll 1\\ \frac{4\Omega^2t}{\gamma} \quad \gamma t \gg 1. \end{cases} \end{align} \tag{ C.6 }$$

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_\gamma = 2\Omega^2/\gamma$ such that $\left\langle x_{q_2} \right\rangle \sim 2D_{\gamma}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² is also given by $\left\langle x^{} \right\rangle_{P_{q^2}} = \sum_j j^2 n(j,t)$, where $n(j,t) = \mathrm{Tr}\{{\rho(t)\pi_j}\}$ can be evaluated numerically (see figure C1 upper panels), or analytically taking the average over the probability distribution associated to π_j

$$\begin{align} P_{\pi_j}\left(x;t\right) = \overline{\delta\left(\left\langle \psi_\xi\left(t\right) | \pi_j | \psi_\xi\left(t\right) \right\rangle-x\right)}. \end{align} \tag{ C.7 }$$

Explicitly we get

$$\begin{align} \left\langle x^{} \right\rangle_{P_{\pi_j}} = n\left(j,t\right) = e^{-\gamma t } \int_{-\pi}^\pi \frac{\mathrm{d} k}{2\pi} e^{-ijk}I_k\left(t\right), \end{align} \tag{ C.8 }$$

notice that, as expected, $\sum_j n(j,t) = 1$, since $I_0(t) = e^{\gamma t}$.