1 Introduction

The time-evolution of closed quantum systems is unitary, deterministic, and governed by Schrödinger equations. By contrast, open quantum systems (see e.g. [1, 2] for reviews) are constantly interacting with their environment. In such cases, dynamics is no longer unitary due to dissipation and mixing effects, and to the flow of information into the (often infinitely many) degrees of freedom of the environment. Within Markovian and weak coupling approximations, such system dynamics are implemented by Lindblad (or Lindblad–Gorini–Kossakowski–Sudarshan) dynamical generators [3, 4], which describe evolution under the assumption that the system-bath interaction is not monitored in any way. The resulting quantum dynamics is deterministic and probability conserving but in general non-unitary.

However, modern experiments can monitor correlations between the system dynamics and the environment through suitable measurement processes [5,6,7,8,9]. For example, a single experiment can yield a time-record of observations, from which the behaviour of the system and the bath can be fully reconstructed. This time-record of events is stochastic because of the fundamental laws of quantum mechanics. It is associated with a quantum trajectory [10,11,12,13] that specifies the evolution of the system state conditioned on the given time-record. Averaging this state over all possible time-records (trajectories) recovers the dynamics generated by a Lindbladian. Going beyond the average, information about dynamical fluctuations is available by analysing stochastic quantum trajectories.

In this paper, we explain how to characterise large dynamical fluctuations of quantum stochastic processes by means of the theory of large deviations (LD) [14,15,16,17,18,19,20,21,22,23,24,25]. In particular, we present results about very general LD functionals which encode information about fluctuations of measurement outcomes. This includes general linear and non-linear functions of the quantum state of the system. We address the two main classes of measurement processes that monitor the interaction of a quantum system with its own environment. One class involves the detection of bath quanta emitted by the system—such as photons or particles—and gives rise to discontinuous quantum jump trajectories [1, 2, 11, 26, 27]. The other involves the continuous monitoring of homodyne currents associated to bath operators, which gives rise to quantum state diffusion [1, 2, 10, 27]. We note that similar equations arise also when describing weak or strong measurements of system observables, see for example [28].

The functionals that we derive and discuss represent the counterpart of level 2.5 LD functionals in jump and diffusion processes [24, 29,30,31,32] for quantum stochastic processes. A short presentation of the results for quantum jump processes has appeared before in [33]. We now present (in Sect. 2) an overview of our main results, including (in Sect. 2.5) an outline of the structure of the following Sections.

2 Outline

2.1 Scope

We consider Markovian open quantum dynamics in which the state of the system is described by a reduced (system) density matrix \(\rho (t)\) that is obtained by tracing out the environment. We focus on finite-dimensional quantum systems described by means of a Hilbert space \(\mathbb {C}^n\), where n is the maximum number of orthogonal (basis) states of the space. The quantum state \(\rho \) is then a Hermitian \(n\times n\) matrix with non-negative eigenvalues and \({\text {Tr}}\rho = 1\). In the Markovian limit, the dynamics of \(\rho \) is given by the Lindblad (or Lindblad–Gorini–Kossakowski–Sudarshan) equation [1,2,3,4]

$$\begin{aligned} {\dot{\rho }} = {{\mathcal {L}}}(\rho ) \end{aligned}$$
(1)

with

$$\begin{aligned} {{\mathcal {L}}}(\rho ) = -i [H,\rho ] + \sum _{i=1}^M\left( L_i \rho L_i^\dag - \frac{1}{2} \left\{ \rho ,L^\dag _i L_i \right\} \right) \, , \end{aligned}$$
(2)

where H is the system Hamiltonian and the \(L_i\) are jump operators that depend on the coupling to the environment. This assumption is sufficiently general to cover the dynamics of several interesting quantum systems in contact with their environment [2]. We sometimes refer to the generator \({{\mathcal {L}}}\) as the Lindbladian.

From the density matrix, it is possible to compute all observable properties of the system. In this work we go further, by considering correlations between system observables and measurements that are made in the environment, as well as time-correlations in the stochastic system dynamics. Two general settings are considered: (i) correlations between system properties and the statistics of quantum jumps, corresponding to emission/absorption (for example of photons) into/from the environment, see Fig. 1a; and (ii) correlations between system properties with the measurement of homodyne currents see Fig. 1b. The theories for these two cases are different in their details, although there are common features. The case of quantum jump detection is discussed in Sect. 3 while homodyne measurements are discussed in Sect. 4.

Fig. 1
figure 1

Sketch of two different experiments for open quantum systems. We show illustrative results for the simple two level system. a Photon counting experiment: a detector reveals the emission (or the absence of emissions) of a photon. The time-record of this measure is a sequence of times at which photons are detected. For each such event detected, the state of the system changes abruptly and collapses to its de-excited state, as it can be seen in the plot of the occupation of the excited state. b Homodyne detection experiment: the output light field emitted by the atom interferes with a local oscillator field and is measured by the detector. The detection outcome consists of the time-integrated homodyne current, measuring the intensity of the quadrature of the light field. The quantum state, in this setting, does not undergo sudden jumps but is instead diffusing as measured by the average occupation

2.2 Unravelling

Some correlations between system and environment can be analysed by tilted variants of Eq. (1), see [34,35,36,37,38,39,40,41,42]. Here we take a different approach, which is to unravel the joint dynamics of the system and the environmental measurements. This enables access to a larger set of dynamical observables and correlations. The theory is based on the stochastic evolution of a pure-state density matrix \(\psi _t\), which is a Hermitian \(n\times n\) matrix for which one eigenvalue is \(+1\) and the others are all zero. This matrix evolves by a stochastic process [2, 10] which we write (schematically) as

$$\begin{aligned} \mathrm {d}\psi _t = b(\psi _t) \mathrm {d}t + \mathrm {d}\omega _t\, , \end{aligned}$$
(3)

where \(\mathrm {d}\omega _t\) represents a random (stochastic) increment for \(\psi \), see below for details. One could also consider the stochastic dynamics of a mixed (non-pure) density matrix. This would be necessary, for instance, in those cases in which the measurement performed on the quantum system can have degenerate outcomes. In these cases, we expect the general theory to be the same, however, some details, such as the space of states in which the stochastic process takes place, would need to be modified (see also Sect. 3.1). Averaging over the noise with suitable initial conditions for \(\psi \), a general result is that

$$\begin{aligned} \mathbb {E}[\psi _t] = \rho (t)\, , \end{aligned}$$
(4)

where \(\mathbb {E}\) is an expectation value for the stochastic process. Hence the (quantum-mechanical) average of any system observable A can be obtained as \(\langle A(t) \rangle = \mathbb {E}[ {\text {Tr}}(A\psi _t)]\).

Note that this construction makes use of a density matrix \(\psi \) that remains normalised at all times \({\text {Tr}}(\psi _t)=1\). Other descriptions of the unravelled dynamics may be expressed in terms of states (or matrices) whose norm (or trace) is also time-dependent [2]. In what follows, we will construct probability distributions for \(\psi _t\), for which it is convenient that this object remains normalised.

Every trajectory of the stochastic process (3) is associated with a time-record for the environmental measurements. For example, if the measurements involve photon counting then the noise \(\omega _t\) causes jumps in \(\psi _t\), and the number of these jumps is, for example, the number of emitted photons. Writing \({N}_t\) for the number of jumps between time 0 and time t, this allows computation of observables such as

$$\begin{aligned} \langle N_t A(\tau ) \rangle = \mathbb {E}[ N_t {\text {Tr}}(A\psi _\tau ) ] \, . \end{aligned}$$
(5)

This is an example of an observable quantity that depends on correlations between the system observable A and the environmental measurement \(N_t\), see for example Ref. [28].

The unravelled system also allows access to objects that are not immediately experimentally observable. In particular, quantum mechanical expectation values are linear functions of \(\psi \) but one may also consider objects that are non-linear. For example, in bipartite systems (decomposed into subsystems A and B), the entanglement of the (pure) density matrix \(\psi \) is

$$\begin{aligned} S_E(\psi ) = -{\text {Tr}}_{\mathrm{A}}( \chi (\psi ) \log \chi (\psi ) )\, , \end{aligned}$$
(6)

where \({\text {Tr}}_{\mathrm{A}}\) denotes a partial trace of subsystem A and \(\chi (\psi ) = {\text {Tr}}_{\mathrm{B}}(\psi )\). Then

$$\begin{aligned} {{\mathcal {S}}}_E(t) = \mathbb {E}\left[ S_E(\psi _t) \right] \end{aligned}$$
(7)

measures the average value of the entanglement shared by the two subsystems. This quantity, obtained as an average over time records, is nowadays receiving a lot of attention [43,44,45,46,47,48,49,50,51].

2.3 Large Deviations at Level 2.5

Our focus in this work is on large deviations of time-integrated quantities. A simple example would be \(N_t\), the number of emitted photons, as above. For large t, the distribution of \(N_t\) is sharply-peaked, in the sense that its mean is proportional to t, while its standard deviation is proportional to \(\sqrt{t}\). Large deviation theory [17,18,19,20,21,22,23,24,25] can be used to analyse the rare events where \(N_t\) differs significantly from its mean value, as \(t\rightarrow \infty \). The statistical properties of these events are described by large deviation theory at level 1, within the classification of Donsker and Varadhan [52,53,54,55].

Here we are concerned with large deviations at a more abstract level of theory, which is called in the LD jargon level 2.5. To explain this, we first consider level 2, which motivates us to define the empirical measure for the (pure-state) density matrix \(\psi \). This is

$$\begin{aligned} \mu _t(\psi ) = \frac{1}{t} \int _0^t \mathrm {d}t' \, \delta ( \psi - \psi _{t'} ) \end{aligned}$$
(8)

where we have introduced a Dirac delta function in the space of density matrices, see later sections for details. For a trajectory of (3), this \(\mu _t(\psi )\) measures (roughly speaking) the fraction of the time interval [0, t] that the system spent in state \(\psi \).

We assume throughout that the system has a unique stationary state, hence for large times \(\mu _t(\psi )\) should converge to some \(P_\infty (\psi )\), which is the steady state distribution for \(\psi \). Large deviation theory at level 2 allows computation of the probability that \(\mu _t\) differs significantly from \(P_\infty \). However, this level 2 theory is not sufficient for our purposes, for example it cannot capture the probability distribution of quantities like \(N_t\). The solution to this problem is to consider the joint statistics of \(\mu _t\) and the empirical fluxes \(Q_t\), which correspond to time-averaged jump rates for all possible quantum jumps. The precise definition of \(Q_t\) depends on the structure of the noise term \(\mathrm {d}\omega _t\) in Eq. (3), see below for details. The level 2.5 theory states that the joint distribution of empirical measure and empirical fluxes behaves as

$$\begin{aligned} \mathrm{Prob}\left[ (\mu _t,Q_t) \approx (\mu ,Q) \right] \asymp \exp [ -t I_{2.5}(\mu ,Q) ] \end{aligned}$$
(9)

where \(I_{2.5}\) is an explicit rate function. The notation here a shorthand which indicates that the random variables \((\mu _t,Q_t)\) should lie inside small sets that contain the values \((\mu ,Q)\), and that the equality is valid on the exponential scale as \(t\rightarrow \infty \). For a rigorous mathematical formulation of LD principles, see for example [16].

Two central results of this paper (following [33]) are explicit formulae for \(I_{2.5}\) for the two classes of Markovian open quantum system that were introduced in Sect. 2.1. These results generalise existing results for large deviations at level 2.5 in classical Markov processes [29,30,31,32, 56, 57].

Large deviation principles (LDPs) at level 2.5 have several applications. Two of the most important are: (i) they give a variational characterisation of large deviations at level-1; and (ii) they allow derivation of general bounds on fluctuations, such as thermodynamic uncertainty relations, see for example [58,59,60,61,62,63,64,65,66]. The connection to level 1 is discussed in detail below; the connection to thermodynamic uncertainty relations was discussed in [33], with a brief recap in Sect. 3.4.3, see also Appendix 4.

2.4 Example Two-State Systems

We illustrate the abstract arguments so far by a simple two-state quantum system, that is \(n=2\). This might represent a single spin or a single qubit. We emphasise that none of our results are restricted to this case, but it is useful for illustrative purposes because it allows a simple representation of the empirical objects \(\mu ,Q\). In this case the most general pure-state density matrix can be represented as

$$\begin{aligned} \psi = \begin{pmatrix} \cos ^2 \frac{\theta }{2}&{} \frac{1}{2}e^{-i\phi } \sin \theta \\ \frac{1}{2}e^{i\phi } \sin \theta &{} \sin ^2\frac{\theta }{2} \end{pmatrix}\, , \end{aligned}$$
(10)

where \((\theta ,\phi )\) are the spherical polar coordinates of a point on the Bloch sphere. This means that the empirical measure \(\mu \) can be interpreted as a probability distribution on this sphere. We briefly describe three types of two-state system, in preparation for the discussion in the rest of the paper.

In a classical two-state system the only possibilities for \(\psi \) correspond to the poles of the Bloch sphere, which are \(\theta =0\) and \(\theta =\pi \). Trajectories for \(\psi \) are restricted to the poles [hence \(b(\psi )=0\) in (3)], and they make discrete jumps from pole to pole, with randomly distributed times. In this case, the empirical measure \(\mu \) always consists of two delta functions at the poles, with weights that indicate the time spent there, see sketch in Fig. 2a. The empirical flux Q is a vector containing two numbers, which are the number of transitions from north to south, and the corresponding number from south to north. The large deviations of \(\mu ,Q\) can be derived from the classical theory for Markov chains at level 2.5 [29,30,31,32, 56, 57].

As a second example we consider a two-state open quantum system where a light source drives transitions between the states, and there is incoherent radiative decay from state \(|\!\!\uparrow \rangle \) to state \(|\!\!\downarrow \rangle \). This corresponds to (2) with \(H=\varOmega \sigma _1\), with \(\sigma _1\) being the first Pauli matrix, \(M=1\), and \(L_1=\sigma ^-\sqrt{\gamma }\) (with \(\sigma ^-=|\!\!\downarrow \rangle \langle \uparrow \!\!|\)) which is also the system pictorially represented in Fig. 1. In this case we explain below that \(\mu \) is supported on a single line on the Bloch sphere, which corresponds to a deterministic evolution given by an effective (non-hermitian) Hamiltonian, starting from the south pole, as in Fig. 2b. The empirical flux Q is a function over the path, parametrised in t by \({\tilde{\psi }}_t\) indicated in Fig. 2b, and provides the rates with which the state of the sytems at the different points in the path has jumped back to the south pole. These jumps correspond to radiative decay events.

Fig. 2
figure 2

Sketches of different empirical measures for a two level system. a In the classical case, where superposition is not possible, the allowed states are only the classical spin configurations with spin pointing up \(|\uparrow \rangle \) (north pole) and with spin pointing down \(|\downarrow \rangle \) (south pole). Therefore, the empirical measure must be given by Dirac deltas in these two points. b In the quantum jump example the state is reset, at every jump, to the south pole and covers a deterministic path \({\tilde{\psi }}_t\) until it jumps again. The empirical measure can thus be solely supported on such a path. c In the quantum state diffusion example, stochastic trajectories are supported on the surface of the Bloch sphere and thus the empirical measure is a function defined over it

Finally, our third example is a two-state open quantum system coupled to a homodyne detector. We consider a fully dissipative dynamics with jump operators \(L_j=i\sigma _j\), with \(j=1,2,3\), proportional to Pauli matrices. The average dynamics is known as the fully dephasing channel, however we discuss how single diffusion trajectories sustain non-zero average coherences at stationarity. In this case Eq. (3) corresponds to diffusion motion of \(\psi \) on the Bloch sphere, and the empirical measure is defined over the whole sphere [see Fig. 2c]. In contrast to the (relatively) simple cases considered so far, the empirical flux Q in this case is a more complicated object: it is related to the empirical current for the spherical diffusion. It turns out, however, that for homodyne quantum trajectories it is also necessary to introduce empirical characterizations of the noises. Details will be discussed below.

2.5 Structure of the Paper

Having set the scene, we outline the structure of what follows. The statistics of quantum jumps are considered in Sect. 3, and those of homodyne currents are discussed in Sect. 4.

These two main Sections are similar in structure: after introductory material in Sects. 3.1 and 3.2 (respectively Sects. 4.1 and 4.2), the level-2.5 LD principles are presented in Sects. 3.3 and 4.3. Then, Sect. 3.4 discusses the relationships between the level-2.5 LD principle and previous results for quantum jumps at level-1, including the quantum Doob transform of [35, 41]. An analogous discussion is given in Sect. 4.4 for homodyne currents. An example system with statistics of homodyne currents is discussed in Sect. 4.5.

Some of the material of Sect. 3 was presented in a shorter form in [33], while that of Sect. 4 is original. Compared to [33], the discussion of Sect. 3 is more comprehensive, particularly in regard to the connections between level-2.5 and level-1, and the similarities and differences between the Doob process of the unravelled system and the quantum Doob transform [35, 41]. The parallel presentation of Sects. 3 and 4 emphasises the general structure of the theory.

3 Quantum Jump Processes

This section discusses the LD properties of quantum stochastic processes where the quantum state makes discontinuous random jumps. For example, such processes can describe experiments where a system emits photons that are detected by some measurement apparatus. When photons are detected, one infers that the system has made a transition into its ground state. A shorter account of these results was presented in Ref. [33], we also review some material from Ref. [35].

3.1 Pure-State Density Matrices and their Calculus

We introduce notation that will be important in the following. Recall that \(\psi _t\) is the pure-state density matrix of the system at time t. This is a Hermitian \(n\times n\) matrix. Denote the set of all Hermitian \(n\times n\) matrices by \({\mathcal {M}}\), and the set of pure-state density matrices by \({{\mathcal {M}}}_{\mathrm{p}}\) (clearly \({{\mathcal {M}}}_{\mathrm{p}}\subset {\mathcal {M}}\)). A generic member of \({{\mathcal {M}}}_{\mathrm{p}}\) has matrix elements

$$\begin{aligned} \psi _{jk} = z^*_j z_k \end{aligned}$$
(11)

where \((z_j)_{j=1}^n\) is a (state) vector with complex elements and \(\sum _{j=1}^n z_j^* z_j=1\). (The notation \(z^*\) indicates the complex conjugate of z.) For stochastic processes evolving mixed-state density matrices, the relevant space of states would be that of positive, unit-trace density matrices, \(\mathcal {M}_{\mathrm{m}}\subset \mathcal {M}\).

The theory that we present is independent of the basis in which the pure state \(\psi \) is represented. However, it is natural to identify a set of classical basis states \(\{ |j\rangle \}_{j=1}^n\) so that \(|j\rangle \) corresponds to a state vector with \(z_j=1\) and \(z_k=0\) for \(k\ne j\). The corresponding matrix \(\psi \) has \(\psi _{jj}=1\) and all other elements are zero, it may be represented as \(\psi = |j\rangle \langle j|\).

Note also that (8) includes a (Dirac) delta function for the matrix \(\psi \). To deal with this we must define integrals over such matrices. It is also useful to define gradients in \({{\mathcal {M}}}\). We achieve this by treating each matrix element as a separate variable, see Appendix 1. For a scalar function \(f=f(\psi )\) [that is, \(f:{{\mathcal {M}}}\rightarrow \mathbb {R}\)] the gradient is a matrix \(\nabla f\) with elements

$$\begin{aligned} (\nabla f)_{jk} = \frac{\partial f}{\partial \psi _{jk}} \; . \end{aligned}$$
(12)

Also, given a matrix \(X\in {{\mathcal {M}}}\) we define

$$\begin{aligned} X\cdot \nabla f = \sum _{jk} X_{jk} (\nabla f)_{jk} \, . \end{aligned}$$
(13)

The theory required for integration is also outlined in Appendix 1. All integrals with respect to \(\psi \) (or \(\psi ')\) are taken over \({{\mathcal {M}}}\) (although in most cases the integrand is non-zero only on \(\mathcal {M}_\mathrm{p}\)). The key results are

$$\begin{aligned} \int d\psi \, \delta (\psi -\chi ) f(\psi ) = f(\chi ) \end{aligned}$$
(14)

and an integration-by-parts formula

$$\begin{aligned} \int d\psi \, X(\psi ) \cdot \nabla f(\psi ) = - \int d\psi f(\psi ) \nabla \cdot X(\psi ) \end{aligned}$$
(15)

where \(X=X(\psi )\) is a matrix-valued function (that is, \(X:{{\mathcal {M}}}\rightarrow {{\mathcal {M}}}\)) and its divergence is \(\nabla \cdot X = \sum _{jk} (\partial X_{jk}/\partial \psi _{jk})\). Since the integrations cover the whole of \({{\mathcal {M}}}\), there are no boundary terms in (15).

3.2 Unravelled Stochastic Dynamics of Open Quantum Systems, and Quantum Jump Trajectories

We now explain how the unravelled quantum dynamics (3) operates in systems with quantum jumps. Recalling Fig. 1, we identify each measurement time-record with a quantum trajectory, which specifies the state of the system at each time t, conditioned on the measurement outcomes obtained. An example of a measurement time-record is the detection plot in Fig. 1a. The time records—and hence the quantum trajectories—are generated by a stochastic process, described by the Belavkin equation [10]

$$\begin{aligned} d{\psi }_t=\mathcal {B}[{\psi }_t]dt+\sum _i\left( \frac{\mathcal {J}_i({\psi }_t)}{{\text {Tr}}[\mathcal {J}_i({\psi }_t)]}-{\psi }_t\right) dn_{it}\, , \end{aligned}$$
(16)

which is the unravelled equation (3), specialised to the jump case. Here, \(d\psi _t\) represents the increment of the pure quantum state \(\psi _t\) in the infinitesimal time interval \([t,t+dt]\) and

$$\begin{aligned} \mathcal {B}[\psi ]=-iH_{\mathrm{eff}}\psi +i\psi H^\dagger _{\mathrm{eff}}-\psi {\text {Tr}}(-iH_{\mathrm{eff}}\psi +i\psi H^\dagger _{\mathrm{eff}})\, , \end{aligned}$$
(17)

with

$$\begin{aligned} H_{\mathrm{eff}}=H-i/2\sum _iL^\dagger _iL_i , \qquad \hbox {and} \qquad \mathcal {J}_i[\psi ]=L_i\psi L^\dagger _i \; . \end{aligned}$$
(18)

The \(dn_{it}\) in (16) are random noise increments whose possible values are 0, 1; they account for the detection events, see Fig. 1a. Only one event can occur in any infinitesimal time period, which means that \(dn_{it}dn_{jt}=\delta _{ij}dn_{it}\); the average noise increment, conditioned on the state being in \(\psi _t\), is \(\mathbb {E}_{\psi _t}[dn_{it}]={\text {Tr}}\left( \mathcal {J}_i[\psi _t]\right) dt\). We emphasise that Eq. (16) describes the time-evolution of a matrix and must be interpreted as a set of equations for increments of matrix elements \((d\psi _{t})_{jk}\).

3.2.1 Comparison of Quantum and Classical Processes

Equation 16 can describe both classical and quantum jump processes, on an equal footing. The relevant classical processes are Markov jump processes over the n classical basis states. They are specified by transition rates W(xy) [from classical state x to the classical state y]. Their trajectories are piecewise-constant: the system remains in a classical state for a (random) time interval before making a discrete jump to some other (classical) state. Hence the allowed values of \(\psi _t\) are the (discrete) classical states \(|j\rangle \langle j|\) with \(1\le j \le n\).

To describe the trajectories of these models one first sets \(H=0\) in (16). Then, for every non-zero rate one introduces a jump operator \( L_{xy}=\sqrt{W(x,y)}|y\rangle \langle x|\). This jump operator generates jumps of \(\psi _t\) from \(|x\rangle \langle x|\) to \(|y\rangle \langle y|\), with rate W(xy). [The indices i in (16) are replaced by indices xy, which label the types of jump.] With these conditions, by starting from a classical configuration state one has a classical state at every later time and \({{\mathcal {B}}}[\psi ]=0\).

Quantum jump processes differ in several important respects from classical processes. First, \(\psi _t\) can include quantum superpositions as well as classical states: this means that \(\psi _t\) can take any value from the set \({{\mathcal {M}}}_{\mathrm{p}}\). Second, trajectories for \(\psi _t\) are piecewise-continuous instead of piecewise constant. In fact, the trajectories are piecewise-deterministic: \(\psi _t\) evolves between jumps as \((\partial \psi _t/\partial t)=\mathcal {B}[{\psi }_t]\) which may be solved as

$$\begin{aligned} \psi _{t+\varDelta t}=\frac{e^{-i \varDelta t H_{\mathrm{eff}}}\, \psi _t \, e^{i \varDelta t H_{\mathrm{eff}}^\dagger }}{{\text {Tr}}\left[ e^{-i \varDelta t H_{\mathrm{eff}}}\, \psi _t \, e^{i \varDelta t H_{\mathrm{eff}}^\dagger } \right] }\, . \end{aligned}$$
(19)

The jumps are discrete, as in the classical case. If a jump occurs at time t via the ith jump operator then the density matrix jumps as

$$\begin{aligned} \psi _t \longrightarrow \frac{\mathcal {J}_i[\psi _t]}{{\text {Tr}}\left( \mathcal {J}_i[\psi _t]\right) }\, . \end{aligned}$$
(20)

This means in particular that while classical jumps occur between discrete configurations, quantum jumps can occur between generic quantum superpositions. Given a system in state \(\psi \), the jump rate into \(\psi '\) (by channel i) is

$$\begin{aligned} w_i(\psi ,\psi ')={\text {Tr}}\left( \mathcal {J}_i[\psi ]\right) \delta \left( \psi ' -\frac{\mathcal {J}_i[\psi ]}{{\text {Tr}}\left( \mathcal {J}_i[\psi ]\right) }\right) \, . \end{aligned}$$
(21)

The \(\delta \) function indicates that the final point of a jump is fully determined by the initial point and the channel.

The fact that the quantum state evolves continuously between jumps also has consequences for the statistics of the times at which the jumps take place. In particular, the probability density function of times between jumps is exponentially distributed in classical jump processes but has a more general structure in quantum systems.

3.2.2 Unravelled Quantum Master Equation

As discussed in Ref. [33], it is useful to derive a dynamical generator that describes the evolution of the quantum state given in (16). (The relevant theory is that of piecewise-deterministic Markov processes [1].) The generator for this process is a linear functional:

$$\begin{aligned} \mathcal {W}[f(\psi )]= \mathcal {B}[\psi ] \cdot \nabla f(\psi ) + \sum _{i}\int d\psi '\, w_i(\psi ,\psi ')\left[ f(\psi ')-f(\psi )\right] \, . \end{aligned}$$
(22)

(If f is a matrix-valued function then \({{\mathcal {W}}}\) acts separately on each matrix element.) The generator has the property

$$\begin{aligned} \frac{d}{dt}\mathbb {E}[ f(\psi _t) ] = \mathbb {E}\left[ {{\mathcal {W}}}[f(\psi _t)] \right] \; . \end{aligned}$$
(23)

We note from (17) that

$$\begin{aligned} {{\mathcal {B}}}[\psi ] = -i[H,\psi ] - \frac{1}{2} \sum _i \left( L^\dag _i L_i \psi + \psi L^\dag _i L_i \right) + \psi \sum _i {\text {Tr}}( L^\dag _i L_i \psi ) \; . \end{aligned}$$
(24)

Hence, taking \(f(\psi )=\psi \) in (22) we find

$$\begin{aligned} \mathcal {W}[ \psi ] = {{\mathcal {L}}}(\psi ) \end{aligned}$$
(25)

where \({{\mathcal {L}}}\) is given by (2). This \({\mathcal {L}}\) is a linear operator. Hence by (23), the time evolution of \(\rho (t)=\mathbb {E}[\psi _t]\) is given by (1).

To avoid any confusion associated with the notation in (25), we discuss briefly the object \({{\mathcal {W}}}[\psi ]\). An alternative notation in (22) would be to write \({{\mathcal {W}}}f\) for the function obtained by operating with \({\mathcal {W}}\) on f, so the left hand side of (22) would be \({{\mathcal {W}}}f(\psi )\). In this case one can define the identity function e by \(e(\psi )=\psi \) and the left hand side of (25) would be \({{\mathcal {W}}}e(\psi )\). Throughout this work, that object is denoted by \({{\mathcal {W}}}[\psi ]\).

Physically, we have shown that averaging the pure state \(\psi _t\) over the trajectories of the unravelled dynamics generates the (mixed) density matrix of the open quantum system of interest. It is a non-trivial feature of these unravelled processes that the expectation value of \(\psi \) obeys a closed equation of motion. (The situation is similar to classical Ornstein–Uhlenbeck processes.)

The process (16) also has a master equation, which is an equation of motion for the probability density for \(\psi _t\), which is denoted by \(P_t(\psi )\). For a generic function f,

$$\begin{aligned} \int d\psi \, f(\psi ) \frac{d}{dt}P_t(\psi ) = \frac{d}{dt}\mathbb {E}[f(\psi _t)] =\int d\psi \, P_t(\psi ) {{\mathcal {W}}}[f(\psi )] \, . \end{aligned}$$
(26)

Since this equation holds for all f, one obtains from (22) that

$$\begin{aligned} \frac{d}{dt}{P}_t(\psi )=-\nabla \cdot \left[ \mathcal {B}[\psi ]P_t(\psi )\right] +\sum _i \int d\psi '\left[ P_t(\psi ')w_i(\psi ',\psi )-P_t(\psi )w_i(\psi ,\psi ') \right] \, , \nonumber \\ \end{aligned}$$
(27)

which is the unravelled quantum master equation [33]. We define an adjoint operator \({{\mathcal {W}}}^\dag \) via \(\int d\psi \, f {{\mathcal {W}}}^\dag [p] = \int d\psi \, p {{\mathcal {W}}}[f]\), which should hold for all pf. Hence from (26) we can also write \((d/dt)P_t(\psi ) = {{\mathcal {W}}}^\dag [P_t(\psi )]\).

Note that (1) is known as the quantum master equation (QME), but the unravelled quantum master equation (27) is a completely different object. In particular, the unravelled QME describes the time-evolution of a probability density function, similar to standard master equations in the theory of stochastic processes. The QME describes the time-evolution of a density matrix, and has a different structure from standard master equations.

3.2.3 Steady State

We assume throughout that the Hamiltonian and jump operators in (16) are such that the process converges for long times to a unique steady state. This means in particular that for any initial condition \(P_0\), the solution of (27) tends to a unique long-time limit which we denote by \(P_\infty \) (see Ref. [67] for conditions on the uniqueness of this invariant measure for quantum Markov chains). The linear operator \({{\mathcal {W}}}\) has eigenvalues which are non-positive, with at least one zero. Since the state space \({{\mathcal {M}}}_{\mathrm{p}}\) is compact, the uniqueness of the steady state means that the zero-eigenvector of \({{\mathcal {W}}}\) is unique and that all other eigenvalues have (strictly) negative real parts. That is, \({{\mathcal {W}}}\) has a positive spectral gap.

The interpretation of \(P_\infty \) is the probability density for \(\psi _t\), in the steady state. We also define the joint probability density \(\varGamma \) for the initial and final points of quantum jumps, in the steady state. This is

$$\begin{aligned} \varGamma _i(\psi ,\psi ') = P_\infty (\psi ) w_i(\psi ,\psi ') \; . \end{aligned}$$
(28)

Also let \(\varGamma \) be a vector whose elements are the \(\varGamma _i\) (for \(1\le i \le M\)).

3.3 LD Principle at Level 2.5

We now formulate the level 2.5 LD principle for these systems, similar to (9). The empirical measure \(\mu _\tau (\psi )\) was defined in (8). It follows from (8,14) that the trajectory-dependent quantity

$$\begin{aligned} \int \mathrm {d}\psi f(\psi ) \mu _\tau (\psi ) = \frac{1}{\tau } \int _0^\tau dt f(\psi _t) \end{aligned}$$
(29)

is the empirical time-average of f. We now define the quantity that plays the role of Q in (9). This is a vector of empirical jump rates, denoted by \(k_\tau \). For a given trajectory, the empirical jump rate for channel i depends on the initial and final points of every jump in the trajectory; it is defined by

$$\begin{aligned} k^i_\tau (\psi ,\psi ')=\frac{1}{\tau }\sum _{\hbox {jumps { j} by channel { i}} }\delta (\psi _j^--\psi )\delta (\psi _j^+ -\psi ') \; \end{aligned}$$
(30)

where the sum is over all the quantum jumps of type (channel) i that occur in the trajectory; the jth jump is from \(\psi _j^-\) to \(\psi _j^+\). Similarly to (29), integrals involving \(k^i_\tau \) generate weighted sums over the jumps: for any function \(g(\psi ,\psi ')\) then

$$\begin{aligned} \int {d}\psi d\psi ' g(\psi ,\psi ') k_\tau ^i(\psi ,\psi ') = \frac{1}{\tau }\sum _{\hbox {jumps { j} by channel { i}} } g(\psi _j^-,\psi _j^+) \; . \end{aligned}$$
(31)

3.3.1 Statement of LD Principle

Since the system has a unique steady state and \({{\mathcal {W}}}\) has a positive spectral gap, it follows that weighted sums of the form (31) converge for large times to fixed (deterministic) values, as do time averages of the form (29). This can be summarised as follows: for \(\tau \rightarrow \infty \) then

$$\begin{aligned} (\mu _\tau ,k_\tau ) \rightarrow (P_\infty ,\varGamma ) \end{aligned}$$
(32)

with probability one (see also [30]).

The LD theory describes rare events where this convergence fails. We state the relevant LD principle before sketching its derivation. The LD principle states that as \(\tau \rightarrow \infty \) then the joint distribution of \((\mu ,k)\) behaves as

$$\begin{aligned} \mathrm{Prob}[\mu ,k]\asymp \exp \left( -\tau I_{2.5}^{\mathrm{qu}}[\mu ,k] \right) \, . \end{aligned}$$
(33)

[This notation has the same meaning as (9), the left hand side is to be interpreted as the probability distribution for \(\mu _\tau ,k_\tau \).]

From (32) one must have \(I_{2.5}^{\mathrm{qu}}[P_\infty ,\varGamma ]=0\). Fixing \((\mu _\tau ,k_\tau )\) specifies the values of all quantities of the form (29,31). This means that the level 2.5 LD principle encodes the (joint) large deviation statistics of all such quantities. The function \(I_{2.5}^{\mathrm{qu}}\) is finite only if the current and flux obey a continuity condition

$$\begin{aligned} \nabla \cdot \left[ \mathcal {B}[\psi ]\mu (\psi )\right] =\sum _i\int d\psi ' \left[ k^i(\psi ',\psi )-k^i(\psi ,\psi ')\right] \, . \end{aligned}$$
(34)

Assuming that this condition holds (and that \(\mu \) is a properly-normalised empirical measure) one has

$$\begin{aligned} I_{2.5}^{\mathrm{qu}}[\mu ,k]=\sum _i \int d\psi d\psi '\, \mathrm{D}\Big [k^i(\psi ,\psi ') \Big | \mu (\psi )w_i(\psi ,\psi ')\Big ] \end{aligned}$$
(35)

where we have introduced the function

$$\begin{aligned} \mathrm{D}[x|y]=x\log (x/y)-x+y\, . \end{aligned}$$
(36)

Equations (3336) fully specify the level 2.5 LD principle for quantum jump trajectories. If the continuity equation (34) does not hold then we set formally \(I_{2.5}^{\mathrm{qu}}[\mu ,k]=+\infty \), this means that \((-1/\tau ) \log \mathrm{Prob}[\mu _\tau ,k_\tau ]\) diverges as \(\tau \rightarrow \infty \).

3.3.2 Derivation of LD Principle

All LD principles in this work are derived by the same general method, based on the Gärtner-Ellis theorem [16, 21]. We first define a moment-generating function (or functional) for the quantity of interest. In this case we consider the empirical measure and flux so we define a generating functional:

$$\begin{aligned} G_\tau [u_1,u_2] = \mathbb {E}\left[ \exp \left( -\tau \int d\psi u_1(\psi ) \mu _\tau (\psi ) - \tau \sum _i \int d\psi d\psi ' \, u_2^i(\psi ,\psi ') k^i_\tau (\psi ,\psi ') \right) \right] \nonumber \\ \end{aligned}$$
(37)

where \(u_1:{{\mathcal {M}}}_{\mathrm{p}}\rightarrow \mathbb {R}\) is a function conjugate to \(\mu \) and similarly \(u_2:{{\mathcal {M}}}_{\mathrm{p}} \times {{\mathcal {M}}}_{\mathrm{p}} \rightarrow \mathbb {R}^M\) is conjugate to k. The corresponding scaled cumulant generating functional (SCGF) is

$$\begin{aligned} \varTheta [u_1,u_2] = \lim _{\tau \rightarrow \infty } \frac{1}{\tau } \log G_\tau [u_1,u_2] \; . \end{aligned}$$
(38)

Then by the Gärtner-Ellis theorem one has (modulo some technical assumptions that are always satisfied in the following):

$$\begin{aligned} I_{2.5}^{\mathrm{qu}}[\mu ,k] = \!\sup _{u_1,u_2} \!\left\{ - \varTheta [u_1,u_2] -\!\int \!d\psi \, u_1(\psi ) \mu (\psi ) - \sum _i \!\int \!d\psi d\psi ' \!u_2^i(\psi ,\psi ') k^i(\psi ,\psi ') \right\} \, .\nonumber \\ \end{aligned}$$
(39)

Moreover, we show in Appendix 2 that \(\varTheta [u_1,u_2]\) may be characterised [33] as the largest eigenvalue of a tilted generator which is a deformed version of \({{\mathcal {W}}}\) in (22):

$$\begin{aligned}&\mathcal {W}_u[f(\psi )]= \mathcal {B}[\psi ] \cdot \nabla f(\psi ) - u_1(\psi ) f(\psi ) \nonumber \\&\quad + \sum _{i}\int d\psi '\, w_i(\psi ,\psi ')\left[ e^{-u_2^i(\psi ,\psi ')} f(\psi ')-f(\psi )\right] \, . \end{aligned}$$
(40)

For many large deviation problems, finding the largest eigenvalue of the tilted generator is prohibitively difficult. However, a key feature of level 2.5 is that the maximisation in (39) can be solved in closed form, yielding (34,35). This computation is described in Appendix 2, it proceeds similarly to that of [30].

3.3.3 Comparison with Level 2.5 for Classical Systems

It is useful to compare the LD principle (33) with corresponding results for classical Markov chains [29, 30], For classical systems as described in Sect. 3.2.1, the empirical jump rate (by channel xy) is simply

$$\begin{aligned} k^{xy}_\tau (\psi ,\psi ')= Q_\tau (x,y) \delta (\psi -|x\rangle \langle x|)\delta (\psi '-|y\rangle \langle y|) \, . \end{aligned}$$
(41)

where \(Q_\tau (x,y)\) is the (classical) empirical jump rate: the number of jumps from the classical state x to the classical state y, normalised by \(\tau \). The corresponding jump rate (21) is

$$\begin{aligned} w_{xy}(\psi ,\psi ')=W(x,y)\delta (\psi -|x\rangle \langle x|)\delta (\psi '-|y\rangle \langle y|)\,. \end{aligned}$$
(42)

Also, the empirical measure \(\mu \) is non-zero only for classical configurations: \(\mu (\psi )=\sum _x\delta (\psi -|x\rangle \langle x|)\mu _{\mathrm{cl}}(x)\) where \(\mu _{\mathrm{cl}}\) is the classical empirical measure, normalised as \(\sum _x \mu _{\mathrm{cl}}(x)=1\). Substituting these facts into \(I_{2.5}^{\mathrm{qu}}(\mu ,k)\) gives

$$\begin{aligned} I_{2.5}^{\mathrm{qu}}(\mu ,k)=\sum _{x\ne y} \left( Q(x,y)\log \frac{Q(x,y)}{\mu _{\mathrm{cl}}(x)W(x,y)}-Q(x,y)+\mu _{\mathrm{cl}}(x)W(x,y)\right) \, , \end{aligned}$$
(43)

which indeed coincides with the classical level 2.5 functional [29, 30]. (The sum runs over pairs of states for which \(W(x,y)\ne 0\).)

To summarise: in the quantum formalism described here, classical jump processes correspond to piecewise constant trajectories for \(\psi _t\), which takes values from a discrete set. In such cases (35) becomes the classical LD principle at level 2.5. The quantum case is more general because \(\psi _t\) follows piecewise-continuous trajectories and can take any value in \({{\mathcal {M}}}_{\mathrm{p}}\).

3.3.4 Auxiliary Process (Doob Transform, Optimally-Controlled Process)

In LD theory, the rate function specifies the probability of rare events. It is also important to characterise the mechanism of these events—that is, the behaviour of trajectories with non-typical values of \((\mu _\tau ,k_\tau )\). The general LD theory explains that these (rare) trajectories can be characterised as typical trajectories of a different system, which we call here the auxiliary process. This Section characterises the auxiliary process associated with the LD result (33).

The derivation is related to a Doob transform and to optimal-control theory, see for example [24, 25]. Note however: the auxiliary process that we describe here is associated to trajectories of the unravelled system, described by a Belavkin equation similar to (16). This is different from the quantum Doob process discussed in [35, 41]. We return to this distinction in Sect. 3.4.3 below.

There is a general recipe for identifying auxiliary processes, using the tilted generator [23]. For any such generator, we define the dominant eigenfunction as the eigenfunction corresponding to the largest eigenvalue. We focus on the tilted generator \({{\mathcal {W}}}_u\), and let \(f_R=f_R(\psi )\) be its dominant eigenfunction. Then the generator of the auxiliary process operates on functions f as

$$\begin{aligned} {{\mathcal {W}}}^A_u[f(\psi )] = f_R(\psi )^{-1} {{\mathcal {W}}}_u[ f(\psi ) f_R(\psi ) ] - \varTheta [u_1,u_2] f(\psi ) \; . \end{aligned}$$
(44)

For \({{\mathcal {W}}}^A_u\) to be a generator of a stochastic process, we require that its largest eigenvalue is zero and that the constant function \(f(\psi )=1\) is the associated eigenvector: \({{\mathcal {W}}}^A_u[1]=0\). This is easily verified for (44). Indeed, this equation allows the auxiliary process to be constructed, dependent on \(u_1,u_2\) and the associated eigenfunction \(f_R\). The generator of the auxiliary process has the same form as (22), but with the rates \(w_i\) replaced by auxiliary rates \(w^A(\psi ,\psi ')\). To find the values of these rates associated to any given \((\mu ,k)\) requires determination of the \(u_1,u_2\) that achieve the maximum in (39). This computation can be performed, formulae for \(w^A\) are given in (147) of Appendix 2. However, the final outcome of the computation can be obtained by direct physical reasoning, as we now explain.

By definition of the auxiliary process, the empirical jump rates k and the empirical measure \(\mu \) are typical of its steady state. This means in particular that the mean jump rate from \(\psi \) to \(\psi '\) must be

$$\begin{aligned} w^A_i(\psi ,\psi ') = \frac{k^i(\psi ,\psi ')}{\mu (\psi )} \; . \end{aligned}$$
(45)

This result fully specifies the auxiliary process for large deviations at level 2.5. It also gives a physical interpretation of the continuity constraint (34): the UQME for the auxiliary process is obtained by replacing w by \(w^A\) in (27). Then (34) says that \(P_t=\mu \) must be a steady state of that equation, consistent with \(\mu \) being the steady state of the auxiliary process.

It is also notable that

$$\begin{aligned} I_{2.5}[\mu ,k] = \sum _i\int d\psi d\psi ' \mu (\psi )\mathrm{D}\Big [w_i^A(\psi ,\psi ')\Big |w_i(\psi ,\psi ')\Big ] \, \end{aligned}$$
(46)

This measures the difference between the auxiliary rates and the original rates of the model. It states that the magnitude of the rate function is determined by the amount by which the rates w must be modified, in order to arrive at a model with the relevant \((\mu ,k)\).

3.4 Full Counting Statistics of Quantum Jumps (LDs at Level-1)

Since the level 2.5 LD principle encodes the probability for large fluctuations of all time-averaged quantities, it can be used to recover the statistics of total quantum jump rates, which are called full counting statistics. We show this explicitly, to indicate how the level 2.5 analysis can be applied. The total (empirical) jump rate for channel i is obtained by integrating the empirical rate \(k^i\) over all initial and final states

$$\begin{aligned} {\bar{k}}^i=\int d\psi d\psi ' k^i(\psi ,\psi ')\, , \end{aligned}$$
(47)

This jump rate obeys a level-1 LD principle, which has been derived in previous work [35, 41] using methods based on tilted Lindblad operators.

This Section shows that the same result can be obtained by contraction from the level-2.5 LD principle, it also explores the relationships between the tilted Lindblad approach and the level-2.5 method described in this work. Specifically, we review the tilted Lindblad method in Sect. 3.4.1, after which Sect. 3.4.2 shows that the same result can be derived from the level 2.5 LD principle. The relationships between the methods are discussed in Sect. 3.4.3, with a focus on the auxiliary process and the quantum Doob process.

3.4.1 Tilted Operator Approach

From (30), the integral (47) is the total number of jumps occurring by channel i in the whole trajectory, normalised by \(\tau \). Also let \({{\bar{k}}}=({\bar{k}}^1,{\bar{k}}^2,\dots ,{\bar{k}}^M)\). For long observation times \(\tau \), the probability distribution of this observable obeys a LD principle

$$\begin{aligned} \mathrm {Prob}({\bar{k}})\asymp \exp \left[ -\tau I_1({\bar{k}}) \right] \, . \end{aligned}$$

To show this, we follow again the general recipe of Sect. 3.3.2. The SCGF is

$$\begin{aligned} \theta _k(\lambda ) = \lim _{\tau \rightarrow \infty } \frac{1}{\tau } \log \mathbb {E}\left[ e^{-\tau \sum _i \lambda _i {\bar{k}}^i_\tau } \right] \end{aligned}$$
(48)

where \({\lambda }=(\lambda _1,\lambda _2,\dots ,\lambda _M)\) is a vector of parameters conjugate to \({\bar{k}}\). The SCGF may be characterised [35] as the largest eigenvalue of a linear operator acting on matrices \(X\in {{\mathcal {M}}}\):

$$\begin{aligned} \mathcal {L}^\dagger _{{\lambda }}(X)=i[H,X]+\sum _i \left( \mathrm{e}^{-\lambda _i} L^\dagger _i X\, L_i -\frac{1}{2}\left[ X L^\dagger _i L_i + L^\dagger _i L_i X \right] \right) \, , \end{aligned}$$
(49)

For \(\lambda =0\) one recovers \({{\mathcal {L}}}^\dag \), which is the adjoint of the operator \({{\mathcal {L}}}\) defined in (2). (This adjoint is defined by the property that \({\text {Tr}}[X{{\mathcal {L}}}(\rho )]={\text {Tr}}[\rho {{\mathcal {L}}}^\dag (X)]\) for all Hermitian matrices \(X,\rho \).) Then \(I_1({\bar{k}})\) can be obtained by Legendre transform [21, 23, 24, 68]

$$\begin{aligned} I_1({\bar{k}})=\sup _{{\lambda }}\left[ -{\bar{k}}\cdot {\lambda }-\theta _k({\lambda })\right] \, . \end{aligned}$$
(50)

(In contrast to level 2.5, neither the SCGF \(\theta _k\) nor the rate function \(I_1\) can be obtained in closed form.)

3.4.2 Level 1 Full-Counting Statistics from the Unravelled Dynamics

We now give a different analysis of full-counting statistics, using the unravelled quantum dynamics (16). The idea is to characterise the SCGF \(\theta _k\) as the largest eigenvalue of a (tilted) generator for the unravelled system, similar to (40). Note that the SCGF \(\theta _k\) in (48) coincides with \(\varTheta [u_1,u_2]\) in (38) if we take \(u_1=0\) and \(u_2^i(\psi ,\psi ')=\lambda _i\). Using (40), it follows that \(\theta _k\) can be characterised as the largest eigenvalue of the tilted generator

$$\begin{aligned} \mathcal {W}_\lambda [f(\psi )]= \mathcal {B}[\psi ] \cdot \nabla f(\psi ) + \sum _{i}\int d\psi '\, w_i(\psi ,\psi ')\left[ e^{-\lambda _i} f(\psi ')-f(\psi )\right] \, . \end{aligned}$$
(51)

We now show explicitly that solving this eigenproblem for \(\theta _k\) is equivalent to finding the largest eigenvalue of (49). To this end, we first show that if \(\theta \) is (any) eigenvalue of \({{\mathcal {L}}}^\dag _\lambda \) then it is also an eigenvalue of \(\mathcal {W}_\lambda \). In this case we have \({{\mathcal {L}}}^\dag _\lambda (\ell ) = \theta \ell \), where \(\ell \) is the relevant eigenmatrix. Now define \(f_\ell (\psi ) = {\text {Tr}}(\ell \psi )\). Since \({{\mathcal {W}}}_\lambda \) is a linear operator we have \({{\mathcal {W}}}_\lambda [ f_\ell (\psi ) ] = {\text {Tr}}(\ell {{\mathcal {W}}}_\lambda [\psi ] )\). Also, it is easily shown [by analogy with (25)] that

$$\begin{aligned} {{\mathcal {W}}}_\lambda [\psi ] = {{\mathcal {L}}}_\lambda (\psi ) \; . \end{aligned}$$
(52)

so that

$$\begin{aligned} {{\mathcal {W}}}_\lambda [ f_\ell (\psi ) ] = {\text {Tr}}\left[ \psi \left( {{\mathcal {L}}}_\lambda ^\dag (\ell ) \right) \right] = \theta f_\ell (\psi ) \end{aligned}$$
(53)

where the second equality uses that \(\ell \) is an eigenmatrix of \({{\mathcal {L}}}^\dag _\lambda \), and the definition of \(f_\ell \). Hence this \(f_\ell \) is an eigenfunction for \({{\mathcal {W}}}_\lambda \) with eigenvalue \(\theta \). However the converse does not hold: there may be eigenvalues of \({{\mathcal {W}}}_\lambda \) that are not eigenvalues of \({{\mathcal {L}}}^\dag _\lambda \).

It therefore remains to show that the largest eigenvalue of \({{\mathcal {W}}}_\lambda \) coincides with the largest eigenvalue of \({{\mathcal {L}}}_\lambda ^\dag \). For a general linear operator, we refer to the eigenfunction corresponding to the largest eigenvalue as the dominant eigenfunction. From our assumption that (16) has a unique steady state, it follows that the dominant eigenfunction of \({{\mathcal {W}}}_\lambda \) is always positive, \(f(\psi )>0\), and that this property is unique to the dominant eigenfunction. Moreover, the theory of Lindblad operators [41] shows that the dominant eigenmatrix \(\ell \) of \({{\mathcal {L}}}_\lambda \) has positive eigenvalues. Since \(\psi \) is a pure state (\(\psi =|z\rangle \langle z|\)) then this implies \(f_\ell (\psi )={\text {Tr}}(\ell \psi )=\langle z|\ell |z\rangle >0\). So \(f_\ell (\psi )\) is an eigenfunction of \({{\mathcal {W}}}_\lambda \) that is always positive – it must be the dominant eigenfunction. Hence the largest eigenvalues of \({{\mathcal {L}}}_\lambda \) and \({{\mathcal {W}}}_\lambda \) are both equal to \(\theta _k\). The level-1 rate function can then be obtained from (50).

Finally, we observe one more way of characterising \(I_1\). By the contraction principle for LDs [16, 21], one has

$$\begin{aligned} I_1({\bar{k}})=\inf _{\mu ,k | {\bar{k}} } I_{2.5}^{\mathrm{qu}}[\mu ,k]\, , \end{aligned}$$
(54)

where the infimum is taken over \((\mu ,k)\), subject to (47). Admissible choices for \(\mu ,k\) in (54) also require that \(\mu \) is normalised and that the continuity condition (34) holds. This minimisation was performed in [33], which verified that it is equivalent to (50). However, the approach here based on the tilted generator \({{\mathcal {W}}}_\lambda \) is a more direct route to the same answer.

3.4.3 Auxiliary Process and Quantum Doob Process

We now turn to the auxiliary process for full-counting statistics, which illustrates the physical connection of the unravelled dynamics to the quantum Doob process of [35, 41], and hence to the tilted Lindblad operator. (The connections are summarized in Fig. 3, below.)

In contrast to the level 2.5 LD principle where explicit results were available, LD results at level-1 rely on the solution to the eigenproblems discussed above. However, the auxiliary rates \(w^A\) are available from (147) [in Appendix 2], in terms of the dominant eigenfunction of \({{\mathcal {W}}}_\lambda \): they are

$$\begin{aligned} w^{\mathrm{A}}_i(\psi ,\psi ')=w_i(\psi ,\psi ')\mathrm{e}^{-\lambda _i}\frac{{\text {Tr}}(\ell \, \psi ')}{{\text {Tr}}(\ell \, \psi )} \;. \end{aligned}$$
(55)

The auxiliary process with these rates reproduces the rare (large deviation) trajectories of the unravelled process, as in Sect. 3.3.4. Similar to (44), the generator of this auxiliary process is

$$\begin{aligned} {{\mathcal {W}}}^A_\lambda [f(\psi )] = {\text {Tr}}(\ell \psi )^{-1} {{\mathcal {W}}}_\lambda [ f(\psi ) {\text {Tr}}(\ell \psi ) ] - \theta _k(\lambda ) f(\psi ) \; . \end{aligned}$$
(56)

In [35, 41], a different kind of auxiliary process was identified, which we call here the quantum Doob process. It corresponds to a Lindblad equation of the form (1), where the Hamiltonian and the jump operators are both modified from the original model of interest. Specifically, the Lindblad generator of this model is given by [35, 41]

$$\begin{aligned} \mathcal {L}_\lambda ^D[\rho ]=\ell ^{1/2} \mathcal {L}_\lambda [\ell ^{-1/2} \rho \ell ^{-1/2}]\ell ^{1/2}-\theta _k(\lambda )\rho \, . \end{aligned}$$
(57)

Using this \(\mathcal {L}_\lambda ^D\) in the Lindblad evolution (1) defines an open quantum system in which the Hamiltonian \({\tilde{H}}\) and the jump operators \({{\tilde{L}}}\) depend on \(\lambda \) as [35, 41]

$$\begin{aligned} {\tilde{H}}&= \frac{1}{2} \ell ^{1/2} \left( H - \frac{i}{2} \sum _i L^\dag _i L_i \right) \ell ^{-1/2} + \hbox {h.c.} \nonumber \\ {\tilde{L}}_i&= \mathrm{e}^{-\lambda _i/2} \ell ^{1/2} L_i \ell ^{-1/2} \end{aligned}$$
(58)

where h.c. denotes the Hermitian conjugate. (We recall that \(\ell \) depends on \(\lambda \).) This new system is the quantum Doob process. It is significant because typical time-records of quantum jumps in the quantum Doob process match exactly the rare time-records that appear as large deviations in the original system. In this sense, the quantum Doob process plays the same role as the auxiliary process for the unravelled dynamics.

Fig. 3
figure 3

Illustration of the relationships between the tilted generators \({{\mathcal {L}}}_\lambda \) and \({{\mathcal {W}}}_\lambda \); the quantum Doob process (described by Lindblad generator \({{\mathcal {L}}}^{D}_\lambda \)); and the unravelled auxiliary processes \({\tilde{\psi }}_t\) and \(\varPsi _t\) (described by classical generators \({{\mathcal {W}}}_\lambda ^{A}\) and \({{\mathcal {W}}}_\lambda ^{D}\)). It is notable that averaging the auxiliary process \({\tilde{\psi }}_t\) does not yield a valid Lindblad evolution for \(\mathbb {E}[{\tilde{\psi }}_t]\). However, the transformation (59) yields an unravelled process \(\varPsi _t\) that is related to the quantum Doob by \(\mathbb {E}[\varPsi _t]=\rho ^{D}\), see Appendix 3

The unravelled dynamics for the quantum Doob process may also be constructed. Let the pure-state density matrix of the auxiliary process of (56) be \({\tilde{\psi }}_t\) and define a new pure-state density matrix:

$$\begin{aligned} \varPsi _t=\frac{\ell ^{1/2} {\tilde{\psi }}_t\ell ^{1/2}}{{\text {Tr}}\left( \ell {\tilde{\psi }}_t\right) }\, . \end{aligned}$$
(59)

The trajectories of this \(\varPsi _t\) define an unravelled jump process which was shown in [33] to coincide with the unravelled dynamics of the quantum Doob process. This is verified in Appendix 3. In other words, the unravelled dynamics of the quantum Doob system can be obtained by deforming the auxiliary process derived here, according to (59). The generator for this unravelled dynamics is denoted by \({{\mathcal {W}}}^D_\lambda \), it can be constructed by analogy with (22), with the transformed Hamiltonian and jump operators from (58) used in place of the original HL.

The relationships between the quantum Doob process and the various unravelled process are illustrated in Fig. 3. It is notable that the auxiliary process described by (55,56) cannot generically be interpreted as the unravelled dynamics of a system obeying Lindblad dynamics (1,2). (The Lindblad form places constraints on the unravelled dynamics which are not satisfied by generic auxiliary processes.) The transformation (59) is essential for relating the unravelled auxiliary processes to the quantum Doob transform.

An application of the level-2.5 LD principle for quantum jumps was considered in [33], which derived a thermodynamic uncertainty relation for photon counts, in the restricted setting of quantum reset processes. An expanded version of that derivation is given in Appendix 4.

3.4.4 Other LDs at Level-1

So far we have considered level-1 LDs of \({{\bar{k}}}\), which are full-counting statistics. These can be investigated either using the unravelled dynamics (via \({{\mathcal {W}}}_\lambda \)) or by a tilted Lindblad operator \({{\mathcal {L}}}_\lambda \). However, working with the unravelled dynamics allows other LD principles at level-1, which cannot be obtained by tilted Lindblad methods.

To see this, consider the fluctuations of a function \(\mathcal {O}=\mathcal {O}(\psi _t)\). Its time-average is

$$\begin{aligned} o_\tau = \frac{1}{\tau } \int _0^\tau dt\, \mathcal {O}(\psi _t)\, , \end{aligned}$$
(60)

In typical cases of interest, the function \(\mathcal {O}(\psi )\) might be the quantum expectation of an operator X, e.g. \(\mathcal {O}(\psi )={\text {Tr}}\left( X\psi \right) \), which is a linear function of \(\psi \). Non-linear functions can also be considered: for example, large deviations of the entanglement entropy of a bipartite quantum system were considered in [33].

The probability density for \(o_\tau \) obeys a LD principle with

$$\begin{aligned} \mathrm{Prob}(o_\tau )\asymp e^{-\tau \phi (o_\tau )}\, , \end{aligned}$$
(61)

where \(\phi (o)\) is the LD rate function. Similar to Sect. 3.4.2, this function may be obtained by contraction from level 2.5. Alternatively the SCGF for \(o_\tau \) is \(\varTheta [u_1,u_2]\) from (38) with \(u_1(\psi )=\lambda {{\mathcal {O}}}(\psi )\) and \(u_2=0\). Hence this SCGF can be obtained as the largest eigenvalue of the appropriate operator \({{\mathcal {W}}}_u\) and the rate function can be obtained by Legendre transform (with \(\lambda \) as the conjugate field). Furthermore, it is also possible to estimate this SCGF by using population dynamics methods [17, 69] applied to the unravelled master equation (27) [70]. We are not aware of any general characterisation of such SCGFs in terms of tilted Lindblad operators.

4 Quantum Diffusion Processes

Many stochastic processes in classical physics are described by differential equations involving Wiener noises (or Langevin equations, or Brownian motions). In the large deviation context, these processes also obey LD principles at level 2.5. The ideas are similar to jump processes, but the technical details are different. In particular the empirical current plays the role of the flux Q in jump processes.

In the quantum context, homodyne measurements on open quantum systems result in random output signals that are related to Brownian motions, recall Fig. 1. (This is in contrast to the photon-detection experiments which are related to jump processes.)

We emphasize that the presentation of this Section is analogous to Sect. 3, with the addition of an example system that is analysed in Sect. 4.5. The LD principle at level-2.5 is presented in Sect. 4.3 and the connection to level 1 is discussed in Sect. 4.4, including the relation between quantum Doob process and unravelled auxiliary process. To set up those results, we briefly review level 2.5 functionals for classical diffusion processes [24, 30], and we explain how these are generalised to the quantum case of homodyne detection experiments [2].

4.1 Summary of LDs at Level 2.5 for Classical Diffusion Processes

As a generic classical diffusion process we take \(x\in \mathbb {R}^d\) evolving by a stochastic differential equation with \(n_\alpha \) independent noises:

$$\begin{aligned} dx_t=A(x_t)dt +\sum _{\alpha =1}^{n_\alpha } B_\alpha dW_t^\alpha \, , \end{aligned}$$
(62)

where \(A:\mathbb {R}^d \rightarrow \mathbb {R}^d\) is a drift term, and \(B_\alpha \in \mathbb {R}^d\) is a vector indicating the strength and direction of noise \(\alpha \) (assumed independent of \(x_t\)). The \(W^\alpha \) are independent Wiener processes with unit variance. We define a diffusion matrix with elements

$$\begin{aligned} D^{ij}_{\mathrm{cl}}=\frac{1}{2} \sum _\alpha B_\alpha ^i B_\alpha ^j \, , \end{aligned}$$
(63)

for \(1\le i,j\le d\). Here and elsewhere, sums over \(\alpha \) are assumed to run from 1 to \(n_\alpha \). The matrix \(D_{\mathrm{cl}}\) is assumed to be invertible [30] which requires (as a necessary condition) that \(n_\alpha \ge d\). We assume that this model has a unique steady state.

In this section, the natural geometry is that of Euclidean space \(\mathbb {R}^d\): gradients such as \(\nabla f\) and dot products such as \(A\cdot \nabla f\) are taken in this space. [In later sections we revert to gradients in \({\mathcal {M}}\), as defined in (12).]

The generator for (62) is \({{\mathcal {W}}}_{\mathrm{diff}}\) which acts on functions \(f:\mathbb {R}^d\rightarrow \mathbb {R}\) as

$$\begin{aligned} {{\mathcal {W}}}_{\mathrm{diff}}[f] =A\cdot \nabla f+\sum _{ij}D_{\mathrm{cl}}^{ij}\frac{\partial ^2 f}{\partial x_i\partial x_j}\, . \end{aligned}$$
(64)

Sums over ij are taken over \(1,2,\dots ,d\). Similar to (26), the generator can be used to derive the Fokker-Planck equation for the time-evolution of the probability density for \(x_t\):

$$\begin{aligned} {\dot{P}}_t=-\nabla \cdot J_{\mathrm{cl}}(P_t) \end{aligned}$$
(65)

where \(J_{\mathrm{cl}}\) is the probability current which depends linearly on \(P_t\). Its elements are

$$\begin{aligned}{}[J_{\mathrm{cl}}(P_t)]_i=A_i P_t-\sum _{j}D_{\mathrm{cl}}^{ij}\frac{\partial P_t}{\partial x_j}\, . \end{aligned}$$
(66)

In the steady state one has \(P_t=P_\infty \) and the associated probability current \(J_{\mathrm{cl},\infty }=J_{\mathrm{cl}}(P_\infty )\) is divergence-free: \(\nabla \cdot J_{\mathrm{cl},\infty }=0\).

4.1.1 Large Deviations

To analyse large deviations at level 2.5, we define the empirical measure using (8), as above. We also define an empirical current as

$$\begin{aligned} {J}^{e}_\tau (x)=\frac{1}{\tau }\int _0^\tau \delta (x_t-x)\circ dx_t \, , \end{aligned}$$
(67)

where the \(\circ \) symbol indicates that the integral uses the Stratonovich convention. For a given trajectory, \(J^e_\tau (x)\) measures the displacement of the system at point x, summed over its visits to that point, and divided by the total time.

For large times we have a result analogous to (32), which is

$$\begin{aligned} (\mu _\tau ,J^e_\tau ) \rightarrow (P_\infty ,J_{\mathrm{cl},\infty }) \end{aligned}$$
(68)

with probability one. The corresponding LD principle is

$$\begin{aligned} \mathrm{Prob} [(\mu _\tau ,J^e_\tau )\approx (\mu ,J)]\asymp \exp \left[ -\tau I_{2.5}^{\mathrm{cl}}(\mu ,J) \right] \, , \end{aligned}$$

which describes the joint statistics of empirical measure and empirical current. The associated rate function is finite only if \(\nabla \cdot J=0\), in which case it takes the value

$$\begin{aligned} I_{2.5}^{\mathrm{cl}}(\mu ,J)=\frac{1}{4}\int dx \left[ J-J_{\mathrm{cl}}(\mu ) \right] \cdot (\mu D_{\mathrm{cl}})^{-1} \left[ J-J_{\mathrm{cl}}(\mu ) \right] \, , \end{aligned}$$
(69)

as discussed in [24, 30].

4.1.2 Empirical Noise

Note the presence of the inverse of \(D_{\mathrm{cl}}\) in (69), so this matrix should be invertible to apply the theory as presented here. For quantum diffusions, the analogue of this matrix may not be invertible. This motivates us to modify the standard theory at level 2.5, as follows. We define an empirical noise

$$\begin{aligned} j^\alpha _\tau (x) =\frac{1}{\tau } \int _0^\tau \delta (x_t-x) dW_t^\alpha \end{aligned}$$
(70)

which is the average noise increment for particles at x. (Note, this integral is taken in the Ito sense, so the empirical noise has mean zero.) Writing j for the vector of empirical noises, it can be shown that the distribution of \((\mu _t,j_t)\) obeys an LDP [similar to (33)]

$$\begin{aligned} \mathrm{Prob} (\mu ,j)\asymp \exp \left[ -\tau I_{\mathrm{noise}}^{\mathrm{cl}}(\mu ,j) \right] \, . \end{aligned}$$
(71)

(We only state this result here, we give the corresponding derivation for the quantum case below. That derivation is easily adapted to this case.) The rate function \(I_{\mathrm{noise}}^{\mathrm{cl}}\) is finite only if a suitable continuity condition \(\nabla \cdot J=0\) holds; this can also be written as

$$\begin{aligned} \nabla \cdot \left[ A \mu + \sum _\alpha B_\alpha j^\alpha - D_{\mathrm{cl}} \nabla \mu \right] = 0 \; . \end{aligned}$$
(72)

The object inside square brackets is the empirical current J for a system with empirical measure \(\mu \) and noise j. It is the sum of \(J_{\mathrm{cl}}(\mu )\) and a term coming from the empirical noise. In cases where the continuity condition (72) holds, the rate function is simply

$$\begin{aligned} I_{\mathrm{noise}}^{\mathrm{cl}}(\mu ,j) = \frac{1}{2} \int dx \sum _\alpha \frac{ j^\alpha (x)^2 }{ \mu (x) } \end{aligned}$$
(73)

The level 2.5 rate function (69) can be obtained from this LD principle by contraction: one minimises \(I^{\mathrm{cl}}_{\mathrm{noise}}\) over all empirical noises \(j^\alpha \) that are consistent with a given empirical current. This minimisation yields the inverse of \(D_{\mathrm{cl}}\) in cases where it exists. (We note in passing that these discussions are not mathematically rigorous, in particular we have not stated precise technical conditions required on (62) in order to obtain this LD principle, although we do insist that the system should have a unique steady state [67]. See also the discussion of the quantum case, below.)

Physically, the meaning of this contraction is that large deviations occur via the least unlikely noise realisations, and \(j_\tau \) characterises these noises. In cases where \(D_{\mathrm{cl}}\) does not have an inverse, there are some empirical currents that cannot be realised by any realisation of the noise. In this case the quantity \(J - J_{\mathrm{cl}}(\mu )\) in (69) is outside the image of \(D_{\mathrm{cl}}\) and the rate function \(I_{2.5}^{\mathrm{cl}}\) is formally infinite.

Note finally, there is a thermodynamic uncertainty principle for currents in these systems [71, 72], it is straightforwardly derived by setting \(\mu =P_\infty \) and \(J=\lambda J_{\mathrm{cl}}(P_\infty )\) in (69), which manifestly solves (72). This construction is valid only if \(D_{\mathrm{cl}}\) has the property that \(J_{\mathrm{cl}} = D_{\mathrm{cl}} F_{\mathrm{cl}}\) may be solved for \(F_{\mathrm{cl}}\) (for all x where \(P_\infty >0\)). The simplest case is when \(D_{\mathrm{cl}}^{-1}\) exists but it is sufficient in general that the steady-state current \(J_{\mathrm{cl}}(P_\infty )\) can be represented as \(\sum _\alpha B_\alpha j^\alpha \), so that there exist realisations of the empirical noise that generate a uniform acceleration of the steady-state current.

4.2 Homodyne Detection Experiments: Unravelled Dynamics

In the case of homodyne detection experiments, quantum systems are monitored by a continuous observation of the quadratures of the bath quantum operators. This type of measurement process allows a detailed characterization of the emissions of the system into the environment [37]. Indeed, since quadrature operators are proportional to the intensity of the light field, the measured homodyne current can provide information not only about the overall number of emitted photons, but also about the nature of the light, i.e. whether this is in a thermal, coherent or more complex state. In ideal conditions, the outcome of a homodyne experiment consists of a record of the time-integrated value of the measured current as a function of time, as shown in Fig. 1b. These time-records are stochastic and depend both on the quantum state and on the specific realization of the noisy interaction between system and environment. In what follows we present a comprehensive discussion of the large deviations in these systems starting from the level 1 statistics for homodyne currents and then deriving the very general level 2.5 functional encoding the statistics of generic observable of the process.

Throughout this section, we consider a system described by the Lindblad evolution (1,2), as for jump processes. However, we slightly change our notation in that jump operators are labelled by m (with \(1\le m \le M\)) instead of by i.

4.2.1 Stochastic Schrödinger Equation

To describe homodyne trajectories, we consider a stochastic Schrödinger equation as in (3). This takes the form of an (Ito) stochastic differential equation similar to (62):

$$\begin{aligned} d \psi _t= \mathcal {L}(\psi _t)dt +\sum _{m=1}^M \mathcal {K}^m(\psi _t) dW_m\, , \end{aligned}$$
(74)

where \({{\mathcal {L}}}\) is the Lindblad operator from (2) and

$$\begin{aligned} \mathcal {K}^m(\psi )&=\kappa _m(\psi )-\psi {\text {Tr}}\left[ \kappa _m(\psi )\right] \, , \nonumber \\ \kappa _m(\psi )&=e^{i\alpha _m}L_m\psi +\psi e^{-i\alpha _m}L_m^\dagger \, , \end{aligned}$$
(75)

where \(\alpha _m\) is a phase factor (see below) and the \(L_m\) are the jump operators appearing in (2). In contrast to (62), the noise strengths \({{\mathcal {K}}}^m\) depend on the state \(\psi \). This means that we must take care to use Ito’s formula when evaluating increments of \(\psi \)-dependent functions. In the literature on quantum diffusions [2], this is implemented by Ito rules

$$\begin{aligned} \mathbb {E}[dW_m]=0, \qquad \mathbb {E}[dW_m dW_n]=\delta _{mn}dt\, . \end{aligned}$$
(76)

The phases \(\alpha _m\) in (75) specify the particular quadrature operator of the environment modes that the experiment is monitoring [37], for each homodyne current. The Lindblad evolution (1) is independent of these phases but the unravelled trajectories can depend qualitatively on the \(\alpha _m\).

Equation (74) is a stochastic differential equation which describes every possible time-record of a homodyne experiment in which the state is being continuously monitored. In particular, a typical outcome consists of the values of the time-integrated homodyne currents \(Q^m\). These are random (trajectory-dependent) quantities, given by

$$\begin{aligned} Q^m_\tau =\int ^\tau _0 d Q^m_t,\quad \text{ with }\quad d Q^m_t={\text {Tr}}\left[ \kappa _m(\psi _t)\right] dt +dW_m\, . \end{aligned}$$
(77)

Let \(Q_\tau \) be a vector whose elements are the \(Q^m_\tau \).

Comparing (74) with (62) one sees that \({\mathcal {L}}\) describes the drift of the diffusion process while \({{\mathcal {K}}}\) describes the noises. From the Ito rules (76) one sees immediately that \(\mathbb {E}(d\psi _t) = \mathbb {E}({{\mathcal {L}}}(\psi _t)) dt\); using that \({{\mathcal {L}}}\) is a linear operator yields \(\mathbb {E}(d\psi _t) = {{\mathcal {L}}}(\mathbb {E}(\psi _t)) dt\). Recalling that \(\mathbb {E}(\psi _t)=\rho _{(t)}\) is the density matrix, one recovers (1). That is, the fact that the drift term is linear in the Ito equation (74) means the expectation value of \(\psi \) obeys a closed (linear) equation. (The same is true for Ornstein–Uhlenbeck equations in the classical setting.)

Unless otherwise stated, we assume in the following that the unravelled process has a unique steady state in which the probability density for \(\psi _t\) is \(P_\infty (\psi )\), as in the case of quantum jump processes.

4.2.2 Unravelled Quantum Master Equation

The next step is to identify the generator for the stochastic process (74). We compute this at the level of the quantum state \(\psi \). Consider a function \(f=f(\psi \)): its increment df in the short time interval \([t,t+dt]\) is obtained by Taylor-expanding to second order:

$$\begin{aligned} d f = \sum _{ij}\frac{\partial f}{\partial \psi _{ij}} (d \psi _t)_{ij}+\frac{1}{2}\sum _{ij,hk}\frac{\partial ^2 f}{\partial \psi _{ij} \partial \psi _{hk}}(d \psi _t)_{ij}(d \psi _t)_{hk} \; \end{aligned}$$
(78)

(It is implicit throughout this section that sums run over all allowed values of the relevant index.) Taking the expectation and using (76) yields

$$\begin{aligned} \mathbb {E}[d f ]= \mathbb {E}\left[ \sum _{ij}\frac{\partial f}{\partial \psi _{ij}} (\mathcal {L}[\psi ])_{ij} +\frac{1}{2}\sum _{ij,hk}\frac{\partial ^2 f}{\partial \psi _{ij}\psi _{hk}}D_{ij,hk}(\psi ) \right] dt \end{aligned}$$
(79)

where

$$\begin{aligned} D_{ij,hk}(\psi )=\sum _m\left( \mathcal {K}^m[\psi ]\right) _{ij} \left( \mathcal {K}^m[\psi ]\right) _{hk}\, \end{aligned}$$
(80)

is the analogue of the classical diffusion matrix \(D_{\mathrm{cl}}\) in this setting (up to a factor of 2). Following that analogy, one sees that that if the number of terms in the sum (M) is not large enough, the matrix D will be degenerate, and the inverse \(D^{-1}\) will not exist. Indeed, this situation is likely to be common for systems under homodyne measurement.

Using (79,23) and recalling (13) we identify the generator for functions of \(\psi \) as

$$\begin{aligned} \mathcal {W}[f(\psi )]= \mathcal {L}[\psi ] \cdot \nabla f(\psi ) +\frac{1}{2}\sum _{ij,hk}D_{ij,hk}(\psi ) \frac{\partial ^2 f}{\partial \psi _{ij}\partial \psi _{hk}}\, , \end{aligned}$$
(81)

which is analogous to the classical result (64). Taking \(f(\psi )=\psi \) recovers again that \(\rho =\mathbb {E}(\psi )\) evolves as in (1).

The analogue of (65) is the unravelled quantum master equation for diffusion processes:

$$\begin{aligned} {\dot{P}}_t=-\nabla \cdot {J}(P_t)\, , \end{aligned}$$
(82)

where \(P_t\) is the probability density for \(\psi \). The corresponding probability current is a matrix-valued function of P, its elements are

$$\begin{aligned} \left[ {J}(P)\right] _{ij}=\mathcal {L}_{ij}P-\frac{1}{2}\sum _{hk} \frac{\partial }{\partial \psi _{hk}}\left( D_{ij,hk} P\right) \; . \end{aligned}$$
(83)

4.3 LD at Level 2.5 for Quantum Diffusions

We now derive a LD principle at level 2.5, following a similar method to Sect. 3.3. The empirical measure is given as usual by (8). Analogous to the classical case from Sect. 4.1, we define the empirical noises

$$\begin{aligned} j^m_\tau (\psi )=\frac{1}{\tau }\int _0^\tau \delta (\psi _t-\psi ) \, dW_m\, . \end{aligned}$$
(84)

Note that

$$\begin{aligned} Q^m_\tau = \tau \int d\psi \left\{ {\text {Tr}}[\kappa ^m(\psi )] \mu (\psi ) + j^m(\psi ) \right\} \; . \end{aligned}$$
(85)

The empirical current is

$$\begin{aligned} {J}^{\mathrm{e}}_\tau (\psi )=\frac{1}{\tau }\int _0^\tau \delta (\psi _t-\psi )\circ d\psi _t\, . \end{aligned}$$
(86)

Note that (86) includes a Stratonovich product, in contrast to the Ito products used elsewhere. Taking care with this fact we show in Appendix 5 that the empirical current is fully determined by the empirical measure and empirical noise, as

$$\begin{aligned}{}[J^e_\tau (\psi )]_{ij} = \mu _\tau (\psi )[{{\mathcal {L}}}(\psi )]_{ij} + \sum _m j_\tau ^m(\psi ) [{{\mathcal {K}}}^m(\psi )]_{ij} - \frac{1}{2} \sum _{hk} \frac{ \partial }{\partial \psi _{hk}} (\mu _\tau (\psi ) D_{ij,hk}(\psi )) .\nonumber \\ \end{aligned}$$
(87)

4.3.1 Large Deviation Principle and Auxiliary Dynamics

We derive a LD principle for \((\mu _\tau ,j_\tau )\) noting that large deviations of \((\mu _\tau ,Q_\tau ,J^e_\tau )\) can then be obtained by contraction. To achieve this, we follow the same steps as Sects. 3.3.2 and 3.3.4. We give a short presentation of the computation, referring to those earlier sections for context and discussion.

Define a moment generating functional for \((\mu ,j)\) that takes as arguments \(a_1:{\mathcal {M}_\mathrm{p}}\rightarrow \mathbb {R}\) and \(a_2:{\mathcal {M}_\mathrm{p}}\rightarrow \mathbb {R}^M\):

$$\begin{aligned} G_\tau [a_1,a_2] = \mathbb {E}\left[ \exp \left( \tau \int d\psi a_1(\psi ) \mu _\tau (\psi ) + \tau \sum _m \int d\psi a_2^m(\psi ) j^m_\tau (\psi ) \right) \right] \, . \end{aligned}$$
(88)

The corresponding SCGF is

$$\begin{aligned} \varTheta [a_1,a_2] = \lim _{\tau \rightarrow \infty } \frac{1}{\tau } \log G_\tau [a_1,a_2] \; . \end{aligned}$$
(89)

The resulting LD principle is

$$\begin{aligned} \mathrm{Prob} (\mu ,j)\asymp \exp \left[ -\tau I^{\mathrm{qu}}_{2.5}(\mu ,j) \right] \, , \end{aligned}$$
(90)

with

$$\begin{aligned} I_{2.5}^{\mathrm{qu}}(\mu ,j) = \sup _{a_1,a_2} \left\{ \int d\psi a_1(\psi ) \mu (\psi ) + \sum _m \int d\psi a_2^m(\psi ) j^m(\psi ) - \varTheta [a_1,a_2] \right\} \; . \end{aligned}$$
(91)

Recall, this last formula should be obtained by applying the Gärtner-Ellis theorem to (89). For a rigorous treatment, this would require technical conditions on \(\varTheta \), which we do not explore here. From a physical perspective, we expect \(\varTheta \) to be well-behaved as long as the unravelled system explores its (unique) steady state within some finite mixing time [67]. We assume that this is the case and the Gärtner-Ellis theorem can be applied – such a requirement is not trivial in systems where D is non-invertible, but we do not expect this to be too restrictive a condition in practice.

This supremum can be computed exactly: we state the (simple) result before outlining the derivation. The rate function is finite only if a continuity equation holds: the empirical current \(J^e\) during large deviation events must converge, for large times, to a current J which must be divergence-free, \(\nabla \cdot {J} = 0\), because the relevant trajectories are stationary. From (87), this requires

$$\begin{aligned} \sum _{ij} \frac{\partial }{\partial \psi _{ij}} \left[ {{\mathcal {L}}}_{ij}\mu + \sum _m {{\mathcal {K}}}^m_{ij} j^m - \frac{1}{2} \sum _{hk} \frac{ \partial }{\partial \psi _{hk}} (\mu D_{ij,hk}) \right] = 0 \; . \end{aligned}$$
(92)

(For compactness of notation, we omit functional dependence on \(\psi \) where this leaves no ambiguity.) In cases where (92) holds then

$$\begin{aligned} {I}^{\mathrm{qu}}_{2.5}(\mu ,j)=\frac{1}{2}\int d\psi \, \sum _{m=1}^M \frac{ j_m^2(\psi ) }{ \mu (\psi ) } \, . \end{aligned}$$
(93)

Just as in the classical case (Sect. 4.1), a thermodynamic uncertainty relation can be derived in this system, if there exist choices of empirical noise such that \( \sum _m {{\mathcal {K}}}^m j^m \) in (92) is proportional to the steady state current \(J(P_\infty )\). One simply substitutes these noises in (93) with \(\mu =P_\infty \) so (92) is easily satisfied. In cases where this construction is not possible, we are not aware of any thermodynamic uncertainty relation.

To derive (92,93), we show in Appendix 5 that \(\varTheta [a_1,a_2]\) from (89) is the largest eigenvalue of the tilted operator

$$\begin{aligned} \mathcal {W}_{a}[f]= \left[ \mathcal {L} +\sum _{m}a_2^m \mathcal {K}^m \right] \cdot \nabla f +\frac{1}{2}\sum _{ij,hk}D_{ij,hk}\frac{\partial ^2 f}{\partial \psi _{ij} \partial \psi _{hk}} + a_1 f +\frac{1}{2}\sum _m (a_2^m)^2 f .\nonumber \\ \end{aligned}$$
(94)

The derivation of (93) from this operator is given in Appendix 5, it is similar to that of Appendix 2 for the jump case.

Similar to Sect. 3.3.4, the auxiliary dynamics is explicit for level 2.5. It may be derived by identifying its generator as

$$\begin{aligned} {{\mathcal {W}}}^A_a[f(\psi )] = f_R(\psi )^{-1} {{\mathcal {W}}}_a[ f(\psi ) f_R(\psi ) ] - \varTheta [a_1,a_2] f(\psi ) \; , \end{aligned}$$
(95)

where \(f_R\) is the dominant eigenvector of \({{\mathcal {W}}}_a\). This is similar to (44). Physically, the meaning of the auxiliary process is that the noise \(dW^m\) develops a (\(\psi \)-dependent) mean value equal to \((j^m/\mu )\). Hence (74) is modified in the auxiliary dynamics as

$$\begin{aligned} d \psi _t= \left[ \mathcal {L}(\psi _t) + \sum _{m} \frac{ \mathcal {K}^m(\psi _t) j^m(\psi _t) }{ \mu (\psi _t) } \right] dt +\sum _{m} \mathcal {K}^m(\psi _t) dW_m\, . \end{aligned}$$
(96)

Similarly, from (77) one sees that the homodyne current in this auxiliary model evolves as

$$\begin{aligned} d Q^m_t= \left[ {\text {Tr}}\left[ \kappa _m(\psi _t)\right] + \frac{j_m(\psi _t)}{\mu (\psi _t)} \right] dt + dW_m \end{aligned}$$
(97)

We recognise (92) as the condition that this auxiliary dynamics has steady-state distribution \(\mu \), similar to the discussion of Sect. 3.3.4 for jump processes.

4.4 LDs of Homodyne Currents (Level 1)

Building on this analysis of level 2.5 LDs, we now consider LDs of homodyne currents. As in the case of jump processes, this will allow us to recover (and extend) earlier work that was based on tilted Lindblad operators [37].

We define the time-averaged homodyne current as

$$\begin{aligned} q^m_\tau = \frac{1}{\tau } Q^m_\tau \; . \end{aligned}$$
(98)

Similar to Sect. 3.4, the level-1 LD principle for this quantity can be computed from the unravelled LD principle at level-2.5, or by using tilted Lindblad operators. The relation between the corresponding auxiliary processes and quantum Doob processes are also analogous to Sect. 3.4. Given that the physical picture is the same as that Section, the presentation here is brief.

4.4.1 Tilted Operators

We analyse large deviations of \(q_\tau \). Its SCGF is

$$\begin{aligned} \theta _q(s) = \lim _{\tau \rightarrow \infty } \frac{1}{\tau } \log \mathbb {E}\left[ \exp \left( - \tau \sum _m s_m q_\tau ^m \right) \right] \end{aligned}$$
(99)

where \(s=(s_1,s_2,\dots ,s_M)\) is the field conjugate to q. By (85), \(\theta _q\) coincides with \(\varTheta [a_1,a_2]\) from (89) with \(a_1=-\sum _m s_m {\text {Tr}}[\kappa ^m(\psi )]\) and \(a_2^m=-s_m\). Hence it suffices to consider the largest eigenvalue of the operator

$$\begin{aligned} \mathcal {W}_{s}[f]&= \left[ \mathcal {L} - \sum _{m} s_m \mathcal {K}^m \right] \cdot \nabla f +\frac{1}{2}\sum _{ij,hk}D_{ij,hk}\frac{\partial ^2 f}{\partial \psi _{ij} \partial \psi _{hk}} \nonumber \\&\quad - \sum _m s_m {\text {Tr}}[\kappa ^m(\psi )] f +\frac{1}{2}\sum _m (s^m)^2 f \, , \end{aligned}$$
(100)

As in Sect. 3.4.2, the dominant eigenfunction f turns out to be linear in \(\psi \). Similar to (25) we have

$$\begin{aligned} {{\mathcal {W}}}_s[\psi ] = {{\mathcal {L}}}_s(\psi ) \end{aligned}$$
(101)

with a tilted Lindblad generator

$$\begin{aligned} \mathcal {L}_s(\rho )=\mathcal {L}(\rho )+\sum _{m} \left[ \frac{1}{2} (s_m)^2 \rho - s_m \kappa _m(\rho ) \right] \; . \end{aligned}$$
(102)

Repeating the argument of Sect. 3.4.2, if \(\ell \) is an eigenmatrix of \( {{\mathcal {L}}}_s^\dag \) with eigenvalue \(\ell \), then \(f(\psi ) = {\text {Tr}}(\ell \psi )\) is an eigenfunction of \({{\mathcal {W}}}_s\), with the same eigenvalue.

The operator \({{\mathcal {L}}}_s\) is a multivariate generalisation of the tilted operator derived in [37]. (The definition of s used here differs from theirs by a factor of 2.) Analysis of this operator shows that the dominant eigenmatrix \(\ell \) has positive eigenvalues, which means that the largest eigenvalue of \({{\mathcal {W}}}_s\) is also the largest eigenvalue of \({{\mathcal {L}}}_s\), and therefore coincides with \(\theta _q(s)\). The operator \({{\mathcal {L}}}_s\) was used in [37] to analyse large deviations of homodyne currents at level-1. Here we have shown how these large deviations can be analysed directly from the unravelled trajectories.

4.4.2 Auxiliary Process and Quantum Doob Process for Homodyne Detection

Similar to the discussion of full-counting statistics in Sect. 3.4.3, one may construct an auxiliary model that reproduces the quantum trajectories associated with large deviation events. One may also construct a quantum Doob-transformed process similar to those in [35, 41].

For the auxiliary process in the unravelled representation, one has from (194) (in Appendix 5) that the empirical noise associated to the rare event is

$$\begin{aligned} j^m(\psi ) = \mu (\psi ) \left( \frac{ {\text {Tr}}[\ell \mathcal{K}^m(\psi )] }{ {\text {Tr}}[\ell \psi ] }- s_m\right) \end{aligned}$$
(103)

(We used that \(a_2^m=-s_m\) and \(f_R={\text {Tr}}[\ell \psi ]\), note that if \(s=0\) then \(\ell \) is the identity and \(j^m=0\), as required.) From (96) the auxiliary process for \(\psi \) is then

$$\begin{aligned} d \psi _t= \left[ \mathcal {L} (\psi _t) + \sum _{m} \mathcal {K}^m(\psi _t) \left( \frac{ {\text {Tr}}[\ell \mathcal{K}^m(\psi _t)] }{ {\text {Tr}}[\ell \psi _t] } - s_m \right) \right] dt + \sum _{m} \mathcal {K}^m(\psi _t) dW_m\, . \end{aligned}$$
(104)

For the quantum Doob process one follows instead the procedure given in [35, 41]. The resulting Lindblad operator is

$$\begin{aligned} \mathcal {L}_s^D[X]=\ell ^{1/2} \mathcal {L}_s[\ell ^{-1/2} X\ell ^{-1/2}]\ell ^{1/2}-\theta (s)X\, . \end{aligned}$$
(105)

This corresponds to a physical model whose typical trajectories allow reconstruction of the homodyne measurement records associated with large deviations event of the original model, analogous to the case of full-counting statistics.

Following again the argument of Appendix 3 [using (101)] shows that the unravelled dynamics of the quantum Doob process can be related to that of the auxiliary process (104) by (59), just as in the case of full counting statistics. (Recall also Fig. 3.)

4.5 Example

We illustrate the application of the level 2.5 formalism with an example from a two-level quantum system, corresponding to a quantum spin-1/2 particle. The space of states in this case admits a pictorial representation in terms of the Bloch sphere. In particular, pure states are parametrised by the spherical polar coordinates \((\theta ,\phi )\) as in (10). We show that the unravelled system corresponds to diffusion on this sphere, and we analyse large deviations in this case. The large deviations of the homodyne currents are quite trivial in this case, so we discuss instead large deviations of the (time-integrated) coherence of the unravelled quantum state.

We consider the dissipative Lindblad dynamics

$$\begin{aligned} \dot{\rho }_t=\sum _{m=1}^3 \left( \sigma _m \rho _t \sigma _m- \rho _t\right) \, , \end{aligned}$$
(106)

where \(\sigma _m\) is the m-th Pauli matrix. The unravelled trajectories are generated by (74), with \(\alpha _m=\pi /2\) in Eq. (75). Hence

$$\begin{aligned} d\psi _t=\sum _{m=1}^3 \left( \sigma _m \psi _t \sigma _m- \psi _t\right) dt +i\sum _{m=1}^3 [\sigma _m,\psi _t ] dW_m \, . \end{aligned}$$
(107)

Note that the stochastic term resembles unitary evolution with a random time-dependent Hamiltonian given (formally) by \(\sum _m \sigma _m (dW_m/dt)\).

4.5.1 Diffusion on the Bloch Sphere

For two-state systems, it is natural to write the density matrix in spherical co-ordinates, as in (10). The equation of motion can then be represented as

$$\begin{aligned} d\begin{pmatrix} \theta _t\\ \phi _t\end{pmatrix}=2\begin{pmatrix} \cot \theta _t\\ 0 \end{pmatrix}dt+\mathcal {K} \begin{pmatrix} dW_1\\ dW_2\\ dW_3 \end{pmatrix} \end{aligned}$$
(108)

with

$$\begin{aligned} \mathcal {K}=2\begin{pmatrix} \sin \phi _t&{}-\cos \phi _t&{}0\\ \cot \theta _t\cos \phi _t&{}\cot \theta _t\sin \phi _t&{}-1 \end{pmatrix}\, . \end{aligned}$$
(109)

This representation emphasises that the set of matrices \({{\mathcal {M}}}_{\mathrm{p}}\) can be parameterised by two real angles. Hence, instead of considering a probability density \(P(\psi )\) as in previous analysis, one may consider a simple probability density for \(\theta ,\phi \). This obeys a Fokker–Planck equation

$$\begin{aligned} \frac{\partial P_t}{\partial t}=-2\frac{\partial }{\partial \theta }\left( \cot \theta P_t\right) +2\frac{\partial ^2P_t}{\partial \theta ^2}+2\csc ^2\theta \frac{\partial ^2P_t}{\partial \phi ^2}\, ; \end{aligned}$$
(110)

alternatively, noting that the uniform distribution on the sphere corresponds to \(P(\theta ,\phi ) = \sin \theta /(4\pi )\) one may define a covariant probability density \(p(\theta ,\phi )=P(\theta ,\phi )/\sin \theta \) which evolves as

$$\begin{aligned} \frac{dp_t}{dt}=2\left[ \frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial p_t}{\partial \theta }\right) +\frac{1}{\sin ^2\theta }\frac{\partial ^2 p_t}{\partial \phi ^2}\right] \, . \end{aligned}$$
(111)

We recognise the right hand side as the Laplacian in spherical co-ordinates. The meaning of this equation is that \(\psi \) undergoes isotropic diffusion on the Bloch sphere, and the steady-state distribution is uniform on the sphere. It follows that the steady-state of the system is time-reversal symmetric, and the probability current in this state is zero.

We work primarily with \(P(\theta ,\phi )\), the probability density for \((\theta ,\phi )\). In this representation, the probability current J is a tangent vector to the Bloch sphere, it has components along the azimuthal (\(\phi )\) and polar (\(\theta \)) directions. That is

$$\begin{aligned} J(P) = ( J_\theta (P), J_\phi (P) ) \end{aligned}$$
(112)

and from (110) one has \(J_\theta = 2[\cot \theta P + (\partial P/\partial \theta )]\) and \(J_\phi = -2 \csc ^2 \theta (\partial P/\partial \phi )\). From this point, the classical large deviation theory of Sect. 4.1 can be applied, as we see below.

We note however that this geometrical representation based on the Bloch sphere is limited to two-state quantum systems (\(n=2\)). The general theory that we have presented is based on probability densities for matrix elements of \(\psi \), which we have denoted by \(P(\psi )\). In this case the probability current is matrix-valued quantity, as in (83). Since such currents may not be intuitive, we give a brief physical discussion of how they appear in this case.

4.5.2 Currents in \({{\mathcal {M}}}\)

The key point from (108) is that while \(\psi \) has 4 elements (two real and two complex), a general increment of \(\psi \) can be described by a two-component vector \((d\theta ,d\phi )\), because \(\psi \) always remains in \({{\mathcal {M}}}_{\mathrm{p}}\). The current is similarly described by the two components \((J_\theta ,J_\phi )\).

There is a useful geometrical structure here: to illustrate it in a simple way, we consider a smooth curve on the Bloch sphere that is described by \(\psi (u)\) with \(0\le u \le 1\). The restriction to smooth paths (and not Brownian motions) avoids complications from Ito’s formula. The curve can be specified in terms of two functions \((\theta (u),\phi (u))\). Then

$$\begin{aligned} \psi '(u) = E_\theta (\psi ) \theta '(u) + E_\phi (\psi ) \phi '(u) \end{aligned}$$
(113)

where primes indicate derivatives, and

$$\begin{aligned} E_\theta (\psi ) = \begin{pmatrix} -\sin \theta &{}\quad \frac{1}{2}e^{-i\phi } \cos \theta \\ \frac{1}{2}e^{i\phi } \cos \theta &{}\quad \sin \theta \end{pmatrix} , \qquad E_\phi (\psi ) = \begin{pmatrix} 0 &{}\quad - \frac{i}{2}e^{-i\phi } \sin \theta \\ \frac{i}{2}e^{i\phi } \sin \theta &{}\quad 0 \end{pmatrix}\nonumber \\ \end{aligned}$$
(114)

are \(\psi \)-dependent matrices which form a basis for the tangent space to the Bloch sphere. It follows that the probability current at the point \(\psi \) is

$$\begin{aligned} J(\psi ) = E_\theta (\psi ) J_\theta (\psi ) + E_\phi (\psi ) J_\phi (\psi ) \end{aligned}$$
(115)

where \(E_\theta ,E_\phi \) are the matrices from (113) and \(J_\theta ,J_\phi \) are the scalar fields defined on the sphere as in (112). The empirical current also has a similar form.

This general picture still holds true for systems with \(n>2\) states: the probability current is a vector field that is everywhere tangent to \(\mathcal {M}_\mathrm{p}\) and can be written as \(J = \sum _\alpha E_\alpha (\psi ) J_\alpha (\psi )\) where the \(E_\alpha \) are matrices that form a basis for the tangent space and the \(J_\alpha \) are real-valued fields. Explicit construction of the basis matrices E and the currents J is not simple for large n (the tangent basis has \(2(n-1)\) independent components). This motivates our general formulation in terms of matrix elements of \(\psi \), which is always applicable.

Fig. 4
figure 4

Large deviation rate function for time-integrated coherence in the diffusion example model (106). The black solid line is the exact LD rate function which has been obtained by Legendre transform of the scaled cumulant generating function, which is the largest real eigenvalue of the twisted Fokker–Planck operator in Eq. (120). The dashed red line is the bound for this rate function which was been found by evalulating the level 2.5 rate function with (122,123)

4.5.3 Large Deviation Analysis

Returning to the example (107), it is easily verified from (77) that \(dQ^m = dW^m\), that is, the three homodyne currents are simple random walks. Hence their rate functions are simple quadratic functions. This is a general feature of systems where the operators \(e^{i\alpha _m}L_m\) are anti-Hermitian. We note that steady state of the Lindblad evolution in our example is \(\rho \propto \mathbf{1}\), the identity matrix, so there are no coherences in the steady state, at this level.

However, while the homodyne currents have Gaussian fluctuations, large deviations of the empirical measure \(\mu \) can have more complex behaviour, and so can coherences within the unravelled system. We consider the coherence of the unravelled density matrix

$$\begin{aligned} {{\mathcal {C}}}(\psi ) = |\psi _{12}| = (\sin \theta )/2 \end{aligned}$$
(116)

This object does have non-trivial fluctuations in the unravelled dynamics. We define its time average

$$\begin{aligned} {{{\bar{c}}}}_\tau =\frac{1}{\tau } \int _0^\tau dt \, \mathcal {C}(\psi _t)\, . \end{aligned}$$
(117)

For large observation times \(\tau \), one has that

$$\begin{aligned} \lim _{\tau \rightarrow \infty } \mathbb {E}[{{{\bar{c}}}}_\tau ]=\int _0^\pi d\theta \int _0^{2\pi } d\phi \, \mathcal {C}(\psi ) P_\infty (\theta ,\phi )=\frac{\pi }{8}\, . \end{aligned}$$
(118)

The coherence also obeys an LD principle \(\mathrm{Prob}({{\bar{c}}}_\tau ) \asymp e^{-\tau I_c({{\bar{c}}}_\tau )}\) with (by a contraction argument)

$$\begin{aligned} I_c({{\bar{c}}}) = \inf _{\mu ,j} I_{2.5}(\mu ,j) \end{aligned}$$
(119)

where the infimum is subject to \(\int d\psi \mu (\psi ) {{\mathcal {C}}}(\psi ) = {\bar{c}}\).

The rate function \(I_c\) can also be obtained as \(\sup _s [ -s{{\bar{c}}} - \varTheta _c(s) ]\) where \(\varTheta _c\) is the largest eigenvalue that solves

$$\begin{aligned} -2\frac{\partial }{\partial \theta }(P \cot \theta )+2\frac{\partial ^2P}{\partial \theta ^2}+2\csc ^2\theta \frac{\partial ^2P}{\partial \phi ^2}+ \frac{s}{2} P \sin \theta = \varTheta _c(s) P \end{aligned}$$
(120)

This problem can be solved numerically (see Fig. 4). We now obtain a simple bound on this rate function via level 2.5. This requires that we construct a suitable \(\mu ,j\) for use in (119).

We use (87) to express the empirical current of our example system in terms of its empirical noises, as

$$\begin{aligned} {J}^e_\theta&=2\mu \cot \theta -2\frac{\partial \mu }{\partial \theta } +2j_1\sin \phi -2j_2\cos \phi \, , \nonumber \\ {J}^e_\phi&=-2\csc ^2\theta \, \frac{\partial \mu }{\partial \phi }+2j_1\cot \theta \cos \phi +2j_2\cot \theta \sin \phi -j_3\, . \end{aligned}$$
(121)

We require \((\partial {J}^e_\theta /\partial \theta ) + (\partial {J}^e_\phi /\partial \phi ) = 0\) in order solve the continuity constraint (92). In fact, the time-reversal symmetry of the original problem and the fact that the coherence is time-reversal symmetric means that the relevant large deviations are realised by trajectories with \({J}^e=0\).

We express \(\mu \) as a probability density for \((\theta ,\phi )\), just like P. We consider a one-parameter family of \((\mu ,j)\), in order to generate a range of values for \({{\bar{c}}}\). Specifically,

$$\begin{aligned} \mu _\lambda (\theta ,\phi )=\frac{1}{2\pi }\frac{e^{-\lambda \sin \theta }\sin \theta }{\int _0^\pi d\theta '\, e^{-\lambda \sin \theta '}\sin \theta '}\, . \end{aligned}$$
(122)

This is independent of \(\phi \), just like \(P_\infty \) (the coherence does not depend on \(\phi \) so there is no reason why the auxiliary process should perturb its distribution away from that of the steady state). For \(\lambda >0\) the distribution \(\mu _\lambda \) is biased towards the poles of the Bloch sphere (less coherence) and for \(\lambda <0\) it biases towards the equator (more coherence). Then suitable empirical noises that achieve \(J^e=0\) are

$$\begin{aligned} j_1= f(\theta ) \sin \phi \, ,\qquad j_2=-f(\theta )\cos \phi \, ,\qquad j_3=0\, , \end{aligned}$$
(123)

with \(f(\theta ) = \frac{\partial \mu }{\partial \theta } -\mu \cot \theta \). Physically, this empirical noise counteracts the tendency of the system to relax towards the steady state \(P_\infty \), leading to trajectories with non-typical values of the coherence. Using these \((\mu _\lambda ,j)\), the value of \(I_{2.5}\) can then be computed. To provide a bound on \(I_c\) one must also compute \(\int d\psi \mu (\psi ) {{\mathcal {C}}}(\psi )\). Performing the integrals numerically, the resulting bound is shown in Fig. 4, together with the (numerically) exact result obtained via (120).

The bound reproduces the general behaviour of \(I_c\) but is not exact. The reason is that the ansatz (122) is not sufficient to fully capture the empirical measure of the rare trajectories that realise the rare event. However, it does have the right general form, especially for small \(\lambda \). The fact that the rare trajectories have \(J^e=0\) is also an accurate reflection of the rare event – one way to see this is that writing (120) as an equation for the covariant density p transforms the eigenvalue problem, into the (time-independent) Schrodinger equation for energy levels of a quantum particle diffusing on a sphere, with potential \((s\sin \theta )/2\). This is a Hermitian eigenvalue problem: these correspond generically to large deviation problems with time-reversal symmetry [19].

5 Outlook

We have analysed unravelled stochastic processes that describe the dynamical evolution of open quantum systems. In particular, we have derived LD principles for these systems at level 2.5, which provide an explicit and general characterisation of the joint fluctuations of the system and environment. We have explained how these LD principles are related to LDs at level 1, which can be analysed by tilted operator methods [35, 37, 42]. We have also discussed the implications of these LD principles for thermodynamic uncertainty relations.

This work opens up a framework for analysing fluctuating behaviour in open quantum systems. It provides a thermodynamic formalism where new interesting nonequilibrium phenomena, such as entanglement phase transitions [43, 44, 46,47,48,49,50] or scrambling of information in stochastic processes [73, 74], can be investigated by means of concepts and tools that proved very powerful in equilibrium statistical mechanics. This makes it possible to look not only at average behaviour but also at the behaviour of higher-order time-correlation functions in quantum stochastic processes.

Since quantum trajectories have a direct relation with measurement outcomes in actual experimental settings involving continuously monitored systems, as for example the photon-counting or homodyne-detection experiments discussed here, this formalism also brings the theoretical investigation of nonequilibrium quantum systems closer to real observations. Given the general connection between large deviation principles and gradient-flow dynamics [75], there are also potential connections of this work to gradient-flow characterisations of Lindblad dynamics [76].