Measuring Fluorescence to Track a Quantum Emitter's State: A Theory Review

TL;DR

This paper reviews the theoretical framework for continuous fluorescence measurement of a qubit, illustrating how different measurement types relate to quantum trajectories and connecting theory with recent experimental results.

Contribution

It provides a unified, Bayesian-based theoretical overview of fluorescence monitoring methods, including photodetection and homodyne detection, and explores advanced topics like the arrow of time in quantum measurement.

Findings

Demonstrates the equivalence of Kraus operator framework and stochastic master equations.

Illustrates methods for quantum trajectory analysis in fluorescence monitoring.

Connects theoretical models with recent experimental observations.

Abstract

We review the continuous monitoring of a qubit through its spontaneous emission, at an introductory level. Contemporary experiments have been able to collect the fluorescence of an artificial atom in a cavity and transmission line, and then make measurements of that emission to obtain diffusive quantum trajectories in the qubit's state. We give a straightforward theoretical overview of such scenarios, using a framework based on Kraus operators derived from a Bayesian update concept; we apply this flexible framework across common types of measurements including photodetection, homodyne, and heterodyne monitoring, and illustrate its equivalence to the stochastic master equation formalism throughout. Special emphasis is given to homodyne (phase-sensitive) monitoring of fluorescence. The examples we develop are used to illustrate basic methods in quantum trajectories, but also to introduce some more advanced topics of contemporary interest, including the arrow of time in quantum measurement, and trajectories following optimal measurement records derived from a variational principle. The derivations we perform lead directly from the development of a simple model to an understanding of recent experimental results.

Equations

$$\rho=\frac{1}{2}\left(\begin{array}[]{cc}1+z&x-iy\\ x+iy&1-z\end{array}\right)=\tfrac{1}{2}\left(\openone+x\hat{\sigma}_{x}+y\hat{\sigma}_{y}+z\hat{\sigma}_{z}\right)$$

$$\begin{split}\left|x+\right\rangle\left\langle x+\right|&=\frac{1}{2}\left(\begin{array}[]{cc}1&1\\ 1&1\end{array}\right)\\ &\neq\frac{1}{2}\left(\begin{array}[]{cc}1&0\\ 0&1\end{array}\right)=\tfrac{1}{2}\left|e\right\rangle\left\langle e\right|+\tfrac{1}{2}\left|g\right\rangle\left\langle g\right|.\end{split}$$

$$\begin{array}{rl}\frac{\hat{h}}{\hslash }=& \underset{\text{qubit}}{\underbrace{{\omega }_{\text{qb}}{\hat{\sigma }}^{+}{\hat{\sigma }}_{-}}}+\underset{\text{single field mode}}{\underbrace{{\omega }_f\left({\hat{a}}^{†}\hat{a}+\frac{1}{2}\right)}}+\underset{\text{interaction }{\hat{h}}_{int}}{\underbrace{\left(g{\hat{\sigma }}^{+}\hat{a}+g^{\ast }{\hat{\sigma }}_{-}{\hat{a}}^{†}\right)}},\end{array}$$

$$\dot{\zeta }=\left(-i{\omega }_{\text{qb}}-\frac{\gamma }{2}\right)\zeta .$$

$${\rho }_{ee}={\zeta }^{\ast }\zeta =(z+1)/2.$$

$${\dot{\rho }}_{ee}=-\gamma {\rho }_{ee}\leftrightarrow \dot{z}=-\gamma \left(z+1\right).$$

$$\dot{\rho }=\underset{\text{Unitary Evol.}}{\underbrace{\frac{i}{\hslash }[\rho ,\hat{H}]}}+\underset{\text{Lindblad Dissipation}}{\underbrace{\sum _c\left({\hat{L}}_c\rho {\hat{L}}_c^{†}-\frac{1}{2}{\hat{L}}_c^{†}{\hat{L}}_c\rho -\frac{1}{2}\rho {\hat{L}}_c^{†}{\hat{L}}_c\right)}},$$

$$\dot{\rho }=\frac{i}{\hslash }[\rho ,\hat{H}]+\gamma {\hat{\sigma }}_{-}\rho {\hat{\sigma }}^{+}-\frac{\gamma }{2}\left({\hat{\sigma }}^{+}{\hat{\sigma }}_{-}\rho +\rho {\hat{\sigma }}^{+}{\hat{\sigma }}_{-}\right),$$

$$\dot{\rho }=\gamma {\hat{\sigma }}_{-}\rho {\hat{\sigma }}^{+}-\frac{\gamma }{2}\left({\hat{\sigma }}^{+}{\hat{\sigma }}_{-}\rho +\rho {\hat{\sigma }}^{+}{\hat{\sigma }}_{-}\right).$$

$$\dot{x}=-\frac{\gamma }{2}x,\dot{y}=-\frac{\gamma }{2}y,\dot{z}=-\gamma (1+z),$$

$$\rho (t+dt)=\frac{{\hat{M}}_r\rho (t){\hat{M}}_r^{†}}{\text{tr}\left({\hat{M}}_r\rho (t){\hat{M}}_r^{†}\right)}.$$

$$|{\psi }_0⟩=\left(\zeta |e⟩+\varphi |g⟩\right)\otimes |0⟩,$$

$$|{\psi }_1⟩=\sqrt{1-ϵ}\zeta |e,0⟩+\varphi |g,0⟩+\sqrt{ϵ}\zeta |g,1⟩.$$

$$\left|\psi_{1}\right\rangle=\left(\begin{array}[]{cc}\sqrt{1-\epsilon}&0\\ \sqrt{\epsilon}\hat{a}^{\dagger}&1\end{array}\right)\left|\psi_{0}\right\rangle\text{ for }\left|\psi_{0}\right\rangle=\left(\begin{array}[]{c}\zeta\\ \phi\end{array}\right)\otimes\left|0\right\rangle,$$

$$\hat{M}_{r}=\left\langle\psi_{r}\right|\left(\begin{array}[]{cc}\sqrt{1-\epsilon}&0\\ \sqrt{\epsilon}\hat{a}^{\dagger}&1\end{array}\right)\left|0\right\rangle,$$

$$\hat{M}_{0}=\left\langle 0\right|\left(\begin{array}[]{cc}\sqrt{1-\epsilon}&0\\ \sqrt{\epsilon}\hat{a}^{\dagger}&1\end{array}\right)\left|0\right\rangle=\left(\begin{array}[]{cc}\sqrt{1-\epsilon}&0\\ 0&1\end{array}\right),$$

$$\hat{M}_{1}=\left\langle 1\right|\left(\begin{array}[]{cc}\sqrt{1-\epsilon}&0\\ \sqrt{\epsilon}\hat{a}^{\dagger}&1\end{array}\right)\left|0\right\rangle=\left(\begin{array}[]{cc}0&0\\ \sqrt{\epsilon}&0\end{array}\right).$$

$$\rho (t+dt)=\frac{{\hat{M}}_1\rho (t){\hat{M}}_1^{†}}{\text{tr}\left({\hat{M}}_1\rho (t){\hat{M}}_1^{†}\right)},$$

$$\rho (t+dt)=\frac{{\hat{M}}_0\rho (t){\hat{M}}_0^{†}}{\text{tr}\left({\hat{M}}_0\rho (t){\hat{M}}_0^{†}\right)}$$

$$\rho (t+dt)=\frac{{\hat{M}}_0\rho (t){\hat{M}}_0^{†}+{\hat{M}}_1\rho (t){\hat{M}}_1^{†}}{\text{tr}\left({\hat{M}}_0\rho (t){\hat{M}}_0^{†}+{\hat{M}}_1\rho (t){\hat{M}}_1^{†}\right)}.$$

$$\dot{\rho }\approx \frac{\rho (t+dt)-\rho (t)}{dt},$$

$$\dot{q}=\mathsf{a}(q)+\mathsf{b}(q)\xi (t);$$

$$d\rho =i[\rho ,\hat{H}]dt+\sum _c\left(\hat{\mathcal{L}}[\rho ,{\hat{L}}_c]dt+\sqrt{{\eta }_c}\hat{\mathcal{M}}[\rho ,{\hat{L}}_c]d{W}_c\right).$$

$$\hat{\mathcal{L}}[\rho ,{\hat{L}}_c]\equiv {\hat{L}}_c\rho {\hat{L}}_c^{†}-\frac{1}{2}\left({\hat{L}}_c^{†}{\hat{L}}_c\rho +\rho {\hat{L}}_c^{†}{\hat{L}}_c\right),$$

$$\hat{\mathcal{M}}[\rho ,{\hat{L}}_c]\equiv {\hat{L}}_c\rho +\rho {\hat{L}}_c^{†}-\rho \text{tr}\left({\hat{L}}_c\rho +\rho {\hat{L}}_c^{†}\right).$$

$$q_{k+1}-q_k=\mathsf{a}({\tilde{q}}_k)\Delta t+\mathsf{b}({\tilde{q}}_k)\Delta W_k,$$

$$d\mathbf{q}=\mathbf{a}(\mathbf{q})dt+\mathbf{b}(\mathbf{q})d\mathbf{W}$$

$$\dot{\mathbf{q}}=\mathbf{A}(\mathbf{q})+\mathbf{b}(\mathbf{q})\mathbf{\xi },$$

$$A_q=a_q-\frac{1}{2}\sum _{j,n}{b}_{jn}{\partial }_n{b}_{jq};$$

$$\dot{\rho }=\hat{\mathcal{L}}[\rho ,{\hat{L}}_X]+\hat{\mathcal{L}}[\rho ,{\hat{L}}_P]+\hat{\mathcal{M}}[\rho ,{\hat{L}}_X]{\xi }_X+\hat{\mathcal{M}}[\rho ,{\hat{L}}_P]{\xi }_P,$$

Full text

I Introduction

The literature on quantum theory and quantum optics is replete with works concerning the spontaneous emission of atoms, across virtually all of its century–long history Einstein (1916); Dirac (1927); Weisskopf and Wigner (1930); Mollow (1969); Ackerhalt et al. (1973); Kimble and Mandel (1976); Milonni (1984); Haroche and Kleppner (1989). The generic case, in which the excited–state population of an emitter decays exponentially on average due to the spontaneous emission of a photon, is a paradigmatic phenomenon in quantum optics. More recently, both the theory Carmichael (1993); Gardiner and Zoller (2004); Percival (1998); Wiseman and Milburn (2010); Barchielli and Gregoratti (2009); Jacobs (2014); Gisin (1984); Barchielli (1986); Caves and Milburn (1987); Diósi (1988); Gisin and Cibils (1992); Gisin and Percival (1992); Wiseman and Milburn (1993a, b, c); Wiseman (1996); Brun (2002); Gough and Sobolev (2004); Jacobs and Steck (2006); Clerk et al. (2010); Santos and Carvalho (2011); Gough et al. (2012); Jordan (2013); Chantasri et al. (2013); Bolund and Mølmer (2014); Chantasri and Jordan (2015); Korotkov (2016); Gross et al. (2018); Gough (2018a, b, 2019, 2020); Crowder et al. (2020) and experiments Gisin et al. (1993); Gambetta et al. (2008); Murch et al. (2013); Weber et al. (2014); Hacohen-Gourgy et al. (2016) about continuous quantum measurement have received considerable attention, and seen rapid progress, revealing new phenomena and insights into the quantum measurement process Hacohen-Gourgy et al. (2016); Lewalle et al. (2017); Dressel et al. (2017); Manikandan and Jordan (2019); Lewalle et al. (2018); Manikandan et al. (2019), and applications to quantum control Patti et al. (2017); Hacohen-Gourgy et al. (2018); Minev et al. (2019). In any such generalized measurement(s) Wiseman and Milburn (2010); Nielsen and Chuang (2000) of some primary system of interest, there must necessarily be some series of interactions between that system and its environment, which allows for information to flow from the primary system to some meter(s) which record the measurement outcome(s) von Neumann (1932). The interaction between the system and environment will necessarily disturb the system of interest in a random way, but inferences about that evolution of the system of interest can be drawn as long as our measurement brings us the information that the environment “learned” by interacting with the system Haroche and Raimond (2006). Generalized measurements can be weak (a small amount of information is acquired about the system state, with correspondingly little disruption to its prior behavior) Aharonov et al. (1988); Tamir and Cohen (2013); Fuchs and Peres (1996); Fuchs and Jacobs (2001); Oreshkov and Brun (2005), or strong (e.g. the system is “collapsed” to an eigenstate of the measurement operator by a projective measurement, such that we have acquired a lot of information at once and disturbed the state by corresponding ramifications in the process) von Neumann (1932); Wiseman and Milburn (2010). Thus, in contrast with the closed–system quantum mechanics described by the Schrödinger equation alone, a physical description of the measurement process necessarily requires that we consider our primary system as being open. A “quantum trajectory” arises when a sequence of measurements are made in time, such that we have a time–series of measurement outcomes, and a corresponding time–series of inferred quantum states of our primary system, based on that information. The process is necessarily stochastic, as there is randomness present in each successive measurement outcome.

Our emphasis will be on tracking the state of a qubit (a two–level quantum system), through its spontaneous emission; this means that the qubit is coupled to a field mode, which is its “environment” in this scheme, and by interrogating this mode in a variety of ways, we will be able to infer a corresponding evolution of the qubit’s state. A quantized electromagnetic field mode is represented by a quantum harmonic oscillator; we will discuss the cases where we interrogate the output field by photodetection (effectively an energy measurement), or by quadrature111The quadrature space of the field is effectively the phase space of the quantum harmonic oscillator describing the field mode in question. In other words, a quadrature is analoguous to the “position” or the “momentum” of a quantum harmonic oscillator, and the product of the noise in orthogonal directions in quadrature space is bounded by the Heisenberg uncertainty principle. A reader unfamiliar with a quadrature phase space representation of a field mode may benefit from perusing e.g. Refs. Leuchs (1988); Silberhorn (2007). measurements (homodyne or heterodyne detection are analoguous to making “position” and/or “momentum” measurements of the oscillator). We are motivated in large part by recent experimental work, in which a superconducting transmon qubit is continuously monitored by homodyne or heterodyne detection of its spontaneous emission, leading to diffusive quantum trajectories Wiseman (1996); Campagne-Ibarcq et al. (2014); Bolund and Mølmer (2014); Jordan et al. (2015); Campagne-Ibarcq et al. (2016a, b); Naghiloo et al. (2016, 2017); Tan et al. (2017); Naghiloo et al. (2020); Ficheux et al. (2018). Following the relevant circuit–QED experiments, the physical setup we have in mind throughout this work involves a single qubit placed inside a cavity, such that microwave photons emitted by the qubit via spontaneous emission are coupled into a transmission line leading to a measurement device. The setup is designed such that photons emitted by the qubit via spontaneous emission are transmitted to the detector, but photons in other modes (e.g. to implement some unitary Rabi rotations on the qubit) are not routed towards it. Such devices allow for high collection efficiency of emitted photons, in contrast with situations in which an atom emits into free space. See Fig. 1 for an illustration, and guides to the experimental details can be found e.g. in Refs. Campagne-Ibarcq et al. (2016a, b); Naghiloo et al. (2016, 2017); Tan et al. (2017); Naghiloo et al. (2020); Ficheux et al. (2018); Ficheux (2018); Naghiloo (2019).

We will proceed by splitting our manuscript into two main parts. In the first part, we describe a qubit open to a decay channel and subsequent measurements from several different perspectives. We carry out the formal treatment of an unmonitored decay channel, using first the typical quantum–mechanical analysis in Sec. II.1, and then by introducing the corresponding master equation in Sec. II.2. We transition towards diffusive trajectories in Sec. III, introducing the Kraus operators that will serve as our primary tool, along with the Stochastic Master Equation (SME). We examine the cases of heterodyne or homodyne detection in detail, in sections IV and V, respectively. We first discuss ideal measurements in which all information is collected, and then describe inefficient measurements in which some information is collected and some is lost in sec. V.2. Each of these steps is represented graphically in Fig. 1. In the second part, we use the examples developed in the first part as a springboard to introduce certain concepts and methods of interest in the current research literature. For example, we are able to use these examples to introduce ideas related to the arrow of time in quantum trajectories in Sec. VI.1, and discuss “optimal paths” (OPs) Chantasri et al. (2013); Chantasri and Jordan (2015); Jordan et al. (2015); Areeya Chantasri (2016); Lewalle et al. (2017, 2018) in Secs. VI.2 and VI.3, which are quantum trajectories which connect given states according to an extremal–probability readout, derived according to a variational principle. Summary, outlook, and further discussion are included in Sec. VII.

II Un–monitored Decay

We review the case of a single qubit whose fluorescence goes unmonitored in two parts. First we review the standard treatment of Weisskopf and Wigner Weisskopf and Wigner (1930); next we introduce an equivalent master equation description of the system Haroche and Raimond (2006); Nielsen and Chuang (2000); Breuer et al. (1997).

A qubit is any two–level quantum system; mathematically–speaking this means that it is described like a spin–$\tfrac{1}{2}$. Physically speaking, a qubit might be any of e.g. a particular transition in an atom or ion, a spin in a quantum dot or diamond nitrogen–vacancy center, or the lowest–two levels of the superconducting “artificial atoms” now used in many experiments. The state of a qubit can be represented as living in the Bloch sphere; we will generically parameterize our single–qubit density matrix with Bloch coordinates according to

[TABLE]

throughout the forthcoming derivations, where $(1+z)/2=\rho_{ee}$ denotes the excited state population. See Fig. 2. It is also necessary that we introduce a distinction between pure states, which can be represented by a state vector $\left|\psi\right\rangle$, and mixed states which require the use of $\rho$. Pure qubit states $\rho=\left|\psi\right\rangle\left\langle\psi\right|$ live on the outer surface of the unit sphere, while more general states $\rho=\sum_{i}\wp_{i}\left|\psi_{i}\right\rangle\left\langle\psi_{i}\right|$ may live inside the sphere; a “mixed” state with more than one non–zero $\wp_{i}$ may be used to describe a system which is imperfectly isolated from its surrounding environment, where the $\wp_{i}$ are effectively probabilities. This description implies that we have a ’classical’ statistical mixture, in which we have a probability $\wp_{i}$ of finding the pure state $\left|\psi_{i}\right\rangle$ in an ensemble; in contrast with a coherent pure state superposition, elements in such a mixture do not interfere with each other. For example, we may consider states which lead to a probability $\tfrac{1}{2}$ for a $\sigma_{z}$–measurement to return $\left|e\right\rangle$ or $\left|g\right\rangle$, such as $\left|x+\right\rangle=\tfrac{1}{\sqrt{2}}(\left|e\right\rangle+\left|g\right\rangle)$. The density matrix for this pure state (which contains the possibility for quantum interference) is not the same as the classical statistical mixture of $\left|e\right\rangle$ and $\left|g\right\rangle$ (where the off–diagonal “coherences” are suppressed), i.e.

[TABLE]

By using the density matrix to describe our state, we may account for all of these options, which is necessary when we have an open system and the possibility of information loss.

Note that $\rho$ is Hermitian ($\rho=\rho^{\dagger}$), and normalization requires that $\text{tr}\left(\rho\right)=1$. We will often represent a qubit’s state and dynamics in terms of the Bloch vector $\mathbf{q}=\{x,y,z\}$ below; such coordinates should be understood in the context of (1). We suppress the “hat” notation on our density operators $\rho$.

II.1 Standard quantum–mechanical treatment

We summarize the typical approach, originally by Weisskopf and Wigner, as a point of departure in describing spontaneous emission. A more complete derivation of the results we summarize can be found in Ref. Jacobs and Steck (2006), and we will stick closely to their conventions for clarity. The Hamiltonian describing the joint system including the qubit and a single mode of the electromagnetic field is of the form

[TABLE]

where $\hbar\omega_{\text{qb}}$ is the energy separation between the two qubit levels of interest, $\omega_{f}$ denotes the frequency of the field mode, and $g$ is a coupling constant between the field and qubit. The first two terms represent the qubit and field, respectively, while the third term describes their interaction. One can regard this model as corresponding to a two level system and quantum harmonic oscillator which are able to exchange excitations. We have raising and lowering operators on the qubit $\hat{\sigma}^{+}=\left|e\right\rangle\left\langle g\right|$ and $\hat{\sigma}_{-}=\left|g\right\rangle\left\langle e\right|$, and on a field mode $\hat{a}\left|n\right\rangle=\sqrt{n}\left|n-1\right\rangle$ and $\hat{a}^{\dagger}\left|n\right\rangle=\sqrt{n+1}\left|n+1\right\rangle$ for Fock states $\left|n\right\rangle$. The interaction term $\hbar g\hat{\sigma}^{+}\hat{a}$ describes the possibility for the atom to become excited by absorbing a photon (which is removed from the field); the adjoint of this term denotes the reverse process, in which the qubit loses energy, emitting a photon which is added to the field mode. We may use the Schrödinger equation to compute the evolution of the joint system and subsequent qubit state amplitudes under the influence of this Hamiltonian. The electro–magnetic environment of the qubit contains many modes, and the apparently incoherent evolution of the qubit associated with spontaneous emission emerges when summing up the action of the couplings to all these modes, each described like $\hat{h}_{int}$ above. Such analysis requires that we average the evolution predicted by the Schrödinger equation over the available density of free–space field modes and sum over polarizations Jacobs and Steck (2006). To good approximation (i.e. the approximations first made by Weisskopf and Wigner), we may simplify the dynamics at timescales much longer than the field periods, eliminating environmental modes from the description. For a generic qubit state $\zeta\left|e\right\rangle+\phi\left|g\right\rangle$, we obtain the evolution of the excited state amplitude

[TABLE]

We have introduced the spontaneous emission rate $\gamma=T_{1}^{-1}$ which is equal to the density of modes coupled resonantly to the qubit via $\gamma=\sum_{j}|g_{j}|^{2}\delta(\omega_{\text{qb}}-\omega_{j})$. Here, $g_{j}$ stands for the strength of the coupling to the $j^{\text{th}}$ mode of the electromagnetic reservoir. The contributions of all the non–resonant modes oscillate and quickly average to zero over a typical time $\tau_{\text{corr}}\ll\gamma^{-1}$; this condition is key in allowing us to obtain a Markovian description of the qubit evolution (4), from which the reservoir dynamics have been completely eliminated. A precise discussion of the approximations leading to Eq. (4) and the order of the associated errors can be found in Ref. Cohen-Tannoudji et al. (1998), Chapter 4.

The excited state population is described by the density matrix element

[TABLE]

It is straightforward to compute that (4) and (5) imply

[TABLE]

The result that spontaneous emission leads to exponential decay of the excited state population at rate $\gamma$ (or with characteristic time $T_{1}=\gamma^{-1}$), absent other dynamics, is among the most fundamental phenomena in the quantum optics literature.

II.2 Master equation treatment

In order to express the complete information about the qubit state at any time in a compact way, and straightforwardly generalize our system (e.g. we might consider coherently driving qubit), it is convenient to formulate the dynamics of spontaneous emission as a master equation for the density operator $\rho$ introduced above. The evolution of an open quantum system in contact with a Markovian environment (i.e. with an environment of very short correlation time $\tau_{\text{corr}}$ with respect to the other time scales in the problem) can, in general, be written as a Lindblad equation; such a master equation is of the form Haroche and Raimond (2006)

[TABLE]

where $\rho$ is the density matrix of the system of primary interest (in this case, the qubit), and each operator $\hat{L}_{c}$ describes a coupling between the system and its environment (in this case, the decay channel). We see a term describing the unitary evolution $\dot{\rho}=\tfrac{i}{\hbar}[\rho,\hat{H}]$, plus the Lindblad term which accounts for information leaking into the environment through any channels to which the system is open.

The arguments of the previous section can be used to show that the case of spontaneous emission corresponds to a single channel characterized by operator $\hat{L}=\sqrt{\gamma}\hat{\sigma}_{-}$ Cohen-Tannoudji et al. (1998), which indicates that the qubit may lose its excitation with an effective coupling rate $\gamma$. The master equation capturing spontaneous emission of a qubit is then

[TABLE]

where the qubit Hamiltonian is $\hat{H}=\hbar\omega_{\text{qb}}\hat{\sigma}_{z}/2$. The unitary part solely induces a rotation at frequency $\omega_{\text{qb}}$ of the qubit state around the $z$-axis of the Bloch sphere, and it is convenient to work in a frame in which this rotation is suppressed. This rotating frame is formally the interaction picture with respect to $\hat{H}$, associated with the transformation $\rho\to e^{i\hat{H}t/\hbar}\rho e^{-i\hat{H}t/\hbar}$. In the following, we will always work in such frame, where the master equation reads

[TABLE]

We can get equations of motion in the Bloch coordinates (in the rotating frame) by computing $\dot{\mathbf{q}}=\text{tr}(\hat{\sigma}_{\mathbf{q}}\dot{\rho})$, yielding

[TABLE]

in perfect agreement with the treatment above (6). The decaying solutions of these equations, initialized from different pure states on the edge of the Bloch sphere, are illustrated in Fig. 3.

For many application, the qubit needs to be driven; to describe such a situation, we must modify the qubit Hamiltonian $\hat{H}$, adding time–dependent terms. In general, the derivation of the master equation describing the dynamics of the qubit’s density operator needs to be carefully redone in presence of this new Hamiltonian, and one may find that the Lindblad term is modified due to the presence of the drive Breuer and Petruccione (2002). Such a modification typically occurs, for example, when the drive causes the qubit to become sensitive to modes of the environment at different frequencies, which are sufficiently separated from each other so as to have a different density of states. For instance, in the case of a quasi–resonant monochromatic drive inducing Rabi oscillations, the emission spectrum of the qubit contains multiple peaks (the famous Mollow triplet Mollow (1969)) separated by the Rabi frequency (which is related to the intensity of the driving). If the environmental density of state varies around the qubit frequency $\omega_{\text{qb}}$ on the Rabi frequency scale, the form of the damping is drastically modified. In Ref. Murch et al. (2012) this effect was exploited to stabilize an arbitrary state of the Bloch sphere. However, provided the drives are weak enough (Rabi frequencies much smaller than the qubit frequency and the inverse correlation time of the reservoir $\tau_{\text{corr}}$), and there is no cavity or other resonance close to the driven qubit’s emission spectrum peaks to cause especially fast variations in the environment spectrum Carmichael (1999), this effect is negligible. Within these conditions, the action of a drive can be simply captured by adding a unitary term $(i/\hbar)[\rho,\hat{H}_{\text{dr}}(t)]$ in Eq. (8); the treatment we develop below assumes this simplest case. While some modification to this simplest scheme may be necessary in adapting it to situations beyond the stated constraints, experiments not explicitly focused on engineering more exotic effects will typically obey these simplifying constraints by default; this simplest scheme we lay out below is thus widely applicable.

III Quantum Trajectories

The treatment of spontaneous emission in the previous section, and in particular the master equation Eq. (9), captures the dynamics of the qubit under the assumption that any information emitted by the qubit (leaking into the environment) during the qubit-field mode interaction is lost forever. We are, however, primarily interested in the case where we, the observer(s), recover some (or ideally all) of this information through measurement(s) on the field mode. In this section, we present the formalism of Kraus operators, which describes the update of the qubit’s state conditioned on acquiring such information.

III.1 Kraus Operator Formalism

The basic idea is that there exist a set of Kraus operators $\hat{M}_{r}$, which describe how the state of our system should be updated, each of them conditioned on acquiring one of the possible measurement outcomes $r$ in the environment during a measurement of duration $dt$, according to Nielsen and Chuang (2000); Wiseman and Milburn (2010); Haroche and Raimond (2006)

[TABLE]

We require that either $\sum_{r}\hat{M}_{r}^{\dagger}\hat{M}_{r}=\openone$, or $\int dr\hat{M}_{r}^{\dagger}\hat{M}_{r}=\openone$, depending on whether the possible measurement outcomes $r$ are discrete or continuous; such a condition tells us that we have a valid (completely positive) transformation on $\rho$, and insures that we have considered a complete, self–consistent set of measurement outcomes. As dictated by the axioms of quantum mechanics, the outcome $r$ is obtained randomly among its possible values based on Born’s rule, which here yields probabilities (or a probability density) $\wp(r|\rho)=\text{tr}\left(\hat{M}_{r}\rho\hat{M}_{r}^{\dagger}\right)$. Note that the denominator of (11), which serves to ensure the updated density matrix is properly normalized, exactly matches this probability. If measurement(s) on the environment is(are) repeated (in our case every $dt$), the successive outcomes and subsequent state updates define a stochastic sequence of states called a quantum trajectory.

Following Ref. Jordan et al. (2015), we construct the particular $\hat{M}$ of interest for the case of a spontaneously–emitting qubit, using a Bayesian probability argument; it is useful to consider a pure state of the qubit and an effective field mode it emits into (initially in vacuum state $\left|0\right\rangle$)

[TABLE]

where $z=2|\zeta|^{2}-1$, and $\left|e\right\rangle$ and $\left|g\right\rangle$ are the excited and ground states of the qubit, respectively. There is a probability $\wp(e)=|\zeta|^{2}$ to find the qubit in $\left|e\right\rangle$ and probability $\wp(g)=|\phi|^{2}$ to find the qubit in $\left|g\right\rangle$, with $\wp(e)+\wp(g)=1$. On phenomenological grounds, we suppose that the probability for an emission event in a time interval $dt$ is given by $\wp(1|e)=\epsilon=\gamma\>dt$, where $\gamma$ is some characteristic rate at which the qubit fluoresces (i.e. $\gamma=1/T_{1}$ is a measurable quantity for a qubit–cavity system). Then $\wp(1|e)\wp(e)=\epsilon|\zeta|^{2}=\wp(e|1)\wp(1)$ and/or $\wp(0|e)\wp(e)=(1-\epsilon)|\zeta|^{2}=\wp(e|0)\wp(0)$ according to Bayes’ theorem. A quantum–coherent state assignment after the short interval $dt$ which reflects these probabilistic considerations is

[TABLE]

In other words, there is some probability for an emission event which involves a photon being created in the output mode ($0\rightarrow 1$), and which shifts qubit population from $\left|e\right\rangle\rightarrow\left|g\right\rangle$, reflecting a common sense understanding of spontaneous emission. Below we will always assume that the measurement time (in practice, a detector integration time) is much faster than the characteristic decay time of the qubit, i.e. we have $dt\ll T_{1}$, or $\epsilon\ll 1$. This is a key condition which will ensure that the quadrature measurements we will eventually consider are weak measurements, and that the subsequent quantum trajectories are diffusive. We also assume that the information we acquire applies to the qubit in real time, which implies that that photon travel time between the qubit and measurement apparatus should be negligible. We may rewrite the change of state from above as

[TABLE]

where $a^{\dagger}$ creates a photon in the relevant cavity/field output ($a^{\dagger}\left|0\right\rangle=\left|1\right\rangle$). The Kraus operators $\hat{M}_{\text{r}}$ in (11) act only on the qubit state, and are obtained by projecting out the field mode in a final state corresponding to some outcome from measuring the field, i.e.

[TABLE]

where $\{\left|\psi_{r}\right\rangle\}$ could be any basis of states of the field mode, which should be chosen based on the kind of measurement being performed and result $r$. All of the examples we consider below rely on a Kraus operator of the form (15). Much of what we do below will revolve around relating different measurements to the appropriate choice of $\left|\psi_{r}\right\rangle$, and then exploring the ramifications that choice has on the measurement backaction and quantum trajectories.

III.2 Photodetection and quantum jump trajectories

As a first example, suppose that we choose our $\left|\psi_{r}\right\rangle$ in the Fock basis, i.e. we consider outcomes of the type $\left|1\right\rangle$ (a photon exits in the field mode in the given timestep), or $\left|0\right\rangle$ (no photon exits), which correspond to making a photodetection measurement. In other words, we imagine counting the photons emitted by the qubit into the field mode, in a time–resolved manner, with a detector integration time $dt\ll\gamma^{-1}$ (equivalently, $\epsilon\ll 1$).

We may define Kraus operators $\hat{M}_{1}$ ($\hat{M}_{0}$) for the single–qubit state update conditioned on a click (no–click) in the detector, according to

[TABLE]

[TABLE]

It is easy to verify that $\hat{M}_{0}^{\dagger}\hat{M}_{0}+\hat{M}_{1}^{\dagger}\hat{M}_{1}=\openone$, such that these measurement operators form a positive operator valued measure (POVM) Nielsen and Chuang (2000). We can say that under continuous photodetection, the qubit state is updated every $dt$ by

[TABLE]

if the detector registers that a photon emerged between $t$ and $t+dt$, or according to

[TABLE]

if no photon reaches the detector. The probability of a click in any given timestep is $\wp_{1}=\text{tr}(\hat{M}_{1}\rho(t)\hat{M}_{1}^{\dagger})=\gamma\>dt\>(1+z)/2$, and the probability of no–click is $\wp_{0}=\text{tr}(\hat{M}_{0}\rho(t)\hat{M}_{0}^{\dagger})$, with $\wp_{0}+\wp_{1}=1$. These expressions reflect the common–sense result that $\wp_{1}$ must vanish when the qubit is in the ground state, i.e. $\wp_{1}=0$ for $z=-1$. Thus, a single quantum trajectory for this photodetection scenario is characterized by a time series of outcomes $r\in\{0,1\}$. Simulation of such a trajectory can be performed by drawing a click/no–click readout from a binomial distribution at each short timestep of duration $dt\ll T_{1}$, and subsequently updating the qubit state $\rho$ according to the appropriate rule above. Results of such a simulation are shown in Fig. 5(a). The trajectories generated by photodetection are an example of “quantum jump” trajectories, for which the qubit state immediately jumps to $\left|g\right\rangle$ when a click event occurs (this is related to the discrete nature of the possible outcomes $r$).

Before moving on to different types of measurements on the output mode, we bridge the gap between our Kraus operator description and the un-monitored decay channel we discussed in the previous section. The situation in which the outcome of the measurement performed on the field mode between $t$ and $t+dt$ is actually unavailable can be captured by averaging the state update over both outcomes, i.e.

[TABLE]

An equation of motion can be obtained by taking

[TABLE]

where the numerator on the RHS is expanded to $O(dt)$. It is then straightforward to verify that the equations (10) reappear exactly, i.e. the procedure just described to obtain $\dot{\rho}$ leads to exactly the same expression as the master equation as described above, and as shown in Fig. 3. A similar procedure allows to show that for any measurement basis $\left|\psi_{r}\right\rangle$ chosen for the field, the master equation is recovered when averaging over all of the outcomes we could have obtained from measurement; we will soon be able to elaborate further on this point.

III.3 Diffusive trajectories and stochastic master equation

In the remainder of this article, we are concerned about measurements on the environment leading to a continuous–valued outcome $r$, e.g. a voltage or current from a detector, leading to “diffusive” trajectories (in contrast with the “jump” trajectories we have just discussed). The specifics of the two most common examples, heterodyne and homodyne measurements, are presented in detail in the following section. Because the evolution during $dt$ is infinitesimal, it is common to write the change in the density operator of the qubit, conditionned on the outcome $r$ obtained at time $t$, under the form of a stochastic master equation (SME); the SME can be seen as an extension of Eq. (7), in which we add a term which accounts for the measurement outcome222Photodetection, as considered above, constitutes a particular “unraveling” of the master equation into stochastic trajectories; the heterodyne and homodyne measurements we subsequently consider are additional possible “unravelings”.. The SME may generically be obtained by expanding an expression of the form (11) to $O(dt)$ Brun (2002); Jacobs and Steck (2006) (detailed examples of this process follow below). The addition of a stochastic element into a differential equation is not trivial, because a genuinely stochastic element is not really differentiable, the way a smooth and well–behaved function is.

Generically, what we will momentarily consider is a type of Langevin equation, or first–order stochastic differential equation of the form

[TABLE]

the term $\mathsf{a}$ is often called the drift term, whereas $\mathsf{b}$ functions as a diffusion constant, and together with the randomly–varying $\xi(t)$, gives stochastic evolution. Equations of this type were first written down to model Brownian motion of small particles Feynman , where complex mechanical forces lead to effectively random kicks in a particle’s position. In our present case, we care about the evolution of a quantum state, and the stochasticity denoted by $\xi(t)$ is a result of the randomness inherent in the quantum measurement process. The particular type of random evolution we consider is delta–correlated Gaussian white noise, obeying $\xi(t)=dW(t)/dt$, where $W(t)$ is called a Wiener process. The Wiener increment $dW(t)=W(t+dt)-W(t)$ is a Gaussian random variable, independent on any past values $dW(s)$ for $s<t$ and characterized by a mean of zero and variance equal to $dt$. These properties lead to a noise term $\xi(t)$ of zero expectation value $\langle\langle\xi(t)\rangle\rangle=0$ and co-variance obeying $\langle\langle\xi(t)\xi(t^{\prime})\rangle\rangle=\delta(t-t^{\prime})$, where the double bracket indicates the ensemble average over realizations of the process. This is suitable for describing the quantum noise arising from measurement in a variety of physical situations, including those we consider below333Strictly speaking, writing $\xi(t)=dW(t)/dt$ is an odd mathematical statement, because $W(t)$ is pure noise and non–differentiable. In practice such substitutions does not cause us a problem in writing down sensible stochastic calculus however. For details, refer e.g. to the books by Gardiner Gardiner and Zoller (2004); Gardiner (2004), or other references on stochastic differential equations, such as Kloeden and Platen (1992).. Some physical justification for the appropriateness of the use of a Gaussian $\xi$ for the examples below is provided in the following sections, and in Appendix A. For $\mathsf{a}=0$ and constant $\mathsf{b}$ (simple diffusion without drift), the variance of an ensemble of diffusing trajectories scales like time; this is summarized by the Itô stochastic calculus rule $dW(t)^{2}=dt$ (or equivalently $\xi(t)^{2}=1/dt$).

The general form of the SME that we use for diffusive quantum trajectories, in units $\hbar\rightarrow 1$, reads Jacobs and Steck (2006); Wiseman and Milburn (2010)

[TABLE]

The super–operators are the Lindblad dissipation term, from (7),

[TABLE]

and the newly–added measurement backaction term

[TABLE]

As before, a Hamiltonian $\hat{H}$ may describe any unitary processes applied to the system (e.g. Rabi drive on a qubit). Each of the operators $\hat{L}_{c}$ describes a particular measurement channel, which is monitored with efficiency $\eta_{c}\in[0,1]$ (where $1$ denotes perfect measurement efficiency, and $\eta_{c}$ is dimensionless). The measurement record associated with any monitored channel, contains one outcome every $dt$, going like $r_{c}\propto\sqrt{\eta_{c}}\langle\hat{L}_{c}+\hat{L}_{c}^{\dagger}\rangle+\xi_{c}$, where the brackets denote the expectation value in state $\rho$. Such an expression for the readout is easy to interpret as a signal $\langle\hat{L}_{c}+\hat{L}_{c}^{\dagger}\rangle$, attenuated due to inefficiency by a factor $\sqrt{\eta_{c}}$, plus quantum noise $\xi_{c}$ intrinsic to the measurement process. A more detailed introductory guide to SME can be found in Ref. Jacobs and Steck (2006).

Channels which are open to the environment, but un-monitored (e.g. typical dephasing mechanisms, or the decay channel in the un-monitored case), can be modeled by placing an operator in the sum over $c$ which is monitored with efficiency $\eta_{c}=0$. The master equation (7) can be recovered from the SME by taking an ensemble average over stochastic trajectories Applying these concepts to the example of a single decay channel introduced above, we see that 1) opening the qubit to an unmonitored decay channel $\hat{L}=\sqrt{\gamma}\hat{\sigma}_{-}$, 2) measuring the qubit fluorescence/decay according to $\hat{L}=\sqrt{\gamma}\hat{\sigma}_{-}$ with efficiency zero, or 3) the average dynamics over an ensemble of stochastic trajectories obtained by continuously monitoring the qubit fluorescence as per $\hat{L}=\sqrt{\gamma}\hat{\sigma}_{-}$, are all equivalent situations. This view from the master equation is also entirely equivalent to that for the Kraus operators, as presented in and around (20). In Fig. 5 we observe this in simulations, observing that the average over many quantum trajectories reproduces the dynamics of the unmonitored case (10), regardless of the character of the individual measurements.

In the sections below, we will formallly compare equations of motion derived from our Kraus operator methods to those from the SME; in order to do this, it is necessary that we briefly comment on a technical issue pertaining to stochastic calculus and the integration of stochastic differential equations. The calculus used to derive and/or manipulate a Langevin of the type above is closely tied to the type of Riemann sum used as the basis of any subsequent integration. If we were integrating an ordinary differential equation, any valid choice of Riemann sum would lead to the same result in the time–continuum limit. This is not so in the stochastic case however; if we suppose that $dW$ is stochastic at every time–scale, different Riemann sums will not converge to the same solutions in the limit anymore! To get the idea, we may consider a discrete update step

[TABLE]

where we have $\tilde{q}_{k}=\beta q_{k+1}+(1-\beta)q_{k}$, and the indices $k$, $k+1$ correspond to times $\Delta t$ apart444Formally, $q(t)-q(0)=\int_{0}^{t}dt^{\prime}\mathsf{a}[q(t^{\prime}),t^{\prime}]+\int_{0}^{t}dW(t^{\prime})\mathsf{b}[q(t^{\prime}),t^{\prime}]$ is more appropriate; the way the integration of the diffusion term, over $dW$, is carried out is both significant and potentially ambiguous. See chapter 4 of Gardiner (2004) for rigorous derivations and more detailed comments.. We highlight two very common conventions: The Itô convention uses $\beta=0$, such that we evaluate drift and diffusion coefficients at the beginning of a timestep, whereas the Stratonovich convention corresponds to $\beta=\tfrac{1}{2}$, such that functions are evaluated according to a trapezoidal rule (see Fig. 4). The form of the SME (23) assumes a derivation based on Itô calculus, in which expansions are made to $O(dt)$ using the rule $dW^{2}=dt$ (i.e. expansions to $O(dt)$ must include explicit expansions to $O(dW^{2})$ in that formalism Jacobs and Steck (2006)). For an accessible and intuitive explanation of this rule, we encourage the interested reader to look at section 4 of Ref. Gough (2018a). Expansions made with regular calculus will lead a Stratonovich equation instead however. In other words, we have to consider two different stochastic calculus conventions, each leading to different differential equations; they give consistent results, however, when paired with the correct integration rules. Specifically: Integrating the equation

[TABLE]

according to the Itô sense ($\beta=0$), is equivalent to performing a Stratonovich integration ($\beta=\tfrac{1}{2}$) on

[TABLE]

where the two drift terms $\mathbf{a}$ and $\mathbf{A}$ are related by the transformation

[TABLE]

$n$ indexes the coordinates (components of $\mathbf{q}$), and $j$ indexes the independent noise(s) on each measurement channel, which are summed. For justification and details see e.g. Gardiner (2004); Kloeden and Platen (1992). We will use this conversion rule to connect different descriptions of the quantum measurement scenarios we consider below.

We make a final remark about numerical simulations before moving on. The appeal of the SME as a theoretical tool is that it expresses quantum trajectory dynamics as a differential equation, similar to how physicists are accustomed to describing classical dynamics; furthermore, the SME readily splits those dynamics into three terms, which make qualitatively distinct contributions to the dynamics. It is worth noting, however, that compared with the case of ordinary differential equations Press et al. (1988), methods for the numerical integration of stochastic differential equations Kloeden and Platen (1992) are more complex, and are accurate only to substantially lower order in $dt$. Additionally, direct numerical integration of the SME does not necessarily preserve the properties of a valid density matrix beyond $O(dt)$, leading to problematic numerical errors unless $dt$ is extremely small; it is consequently numerically preferable to execute simulations of stochastic quantum trajectories by direct application of a positive mapping, as in (11) or similar, when possible. The interested reader may find further comments in this vein e.g. in Ref. Rouchon and Ralph (2015).

IV Single–Qubit Heterodyne Trajectories

We now begin looking at diffusive quantum trajectories due to heterodyne detection. What follows is essentially a review of the simplest non–trivial case described more extensively in Ref. Jordan et al. (2015), and corresponding to the experimental implementation e.g. of Ref. Campagne-Ibarcq et al. (2016a). In the language of quantum–limited amplifiers (QLAs), which are essential to realizing experiments involving individual quantum trajectories, our meaning of “heterodyne” corresponds to “phase–preserving” amplification (e.g. see Clerk et al. (2010); Bergeal et al. (2010); Emmanuel Flurin (2014); Roy and Devoret (2016) or similar, regarding implementations in circuit QED scenarios). See Fig. 1. Owing to the mixing of the fluorescence signal with a coherent state of the field (the “local oscillator”, or LO), the heterodyne measurement gives access to both quadratures of the field, with a symmetric uncertainty. A reader unfamiliar with a quadrature phase space representation of a field mode may benefit from perusing e.g. Ref. Silberhorn (2007). When performed with an ideal QLA, this scheme is formally equivalent to projecting the field mode into the basis of the coherent states Jordan et al. (2015).

IV.1 Stochastic Master Equation Treatment

The SME is given in Eq. (23), and provides one of the most–used approaches to modeling diffusive quantum trajectories arising from continuous weak measurement Brun (2002); Jacobs and Steck (2006). We will consider an idealized measurement in the rotating frame, characterized by $\hat{H}=0$ (no unitary dynamics), $\hat{L}_{X}=\hat{\sigma}_{-}\sqrt{\gamma/2}$, and $\hat{L}_{P}=i\hat{\sigma}_{-}\sqrt{\gamma/2}$, where there is no dephasing channel and the measurement efficiency $\eta=1$ is perfect. We can make qualitative sense of the two operators $\hat{L}_{X}$ and $\hat{L}_{P}$ by understanding that $\hat{\sigma}_{-}$ indicates that our measurement is being made through a decay channel, and that $\hat{L}_{X}$ and $\hat{L}_{P}$ are associated, respectively, with the information encoded in the two quadratures $\hat{X}$ and $\hat{P}$ of the field read out by the heterodyne measurement; the factor $i$ between $\hat{L}_{X}$ and $\hat{L}_{P}$ is the $90^{\circ}$ phase between these two orthogonal directions in the $XP$–plane (often also conventially labeled as the $IQ$–plane).

The resulting SME is then

[TABLE]

where $\hat{\mathcal{L}}$ and $\hat{\mathcal{M}}$ are still the Lindblad dissipation, and measurement backaction terms, respectively. The Gaussian white noise for the measurement channels is characterized by each $\xi(t)\sim dW/dt$. We may obtain equations of motion in terms of Bloch sphere coordinates using $\dot{\mathbf{q}}=\text{tr}(\hat{\sigma}_{\mathbf{q}}\>\dot{\rho})$, yielding

[TABLE]

in agreement with the result in eq. (25) of Jordan et al. (2015) (for $u=1+z$, $\eta=1$, and $\gamma_{\phi}=0$, in their notation). The stochastic readouts (signals arising from the measurement process) are given by

[TABLE]

Notice that the average path given by these equations (where averages over an ensemble lead to $\xi\rightarrow 0$, since these are zero–mean stochastic variables) obeys the same basic fluorescence relations (10). This is a typical example of the relationship between an un-monitored and continuously–monitored system, as we have discussed in general above.

We will interpret the equations (31) as being equations suitable for Itô integration and stochastic calculus (consistent with the assumptions used to derive (23) in the first place Jacobs and Steck (2006)). It will also be useful to have the corresponding Stratonovich version of this system of equations, which can be manipulated using regular calculus. In this case, the conversion (29) can be written as

[TABLE]

a trio of Stratonovich equations corresponding to the Itô equations (31) are obtained by substituting this new drift vector (33) into (28).

IV.2 Kraus Operator Treatment

We now consider the corresponding Kraus operator treatment of this situation. As discussed previously, a heterodyne measurement effectively projects the fluorescence signal onto a coherent state (77) at each measurement timestep, such that we write down an operator

[TABLE]

We will use a substitution for the readouts given by

[TABLE]

the prefactor $\sqrt{dt/2}$ is chosen because it generates statistics consistent with the shot–noise of the coherent state LO; for clarification see Jordan et al. (2015) and/or appendix A. With this substitution, we have a Kraus operator

[TABLE]

which may be used to update the state (using (11) with $\hat{M}_{r}\rightarrow\hat{M}_{\alpha}$) conditioned on acquiring a measurement record drawn from the probability density $\wp(r_{X},r_{P}|\rho(t))=\mathcal{N}\text{tr}\left(\hat{M}_{\alpha}\rho(t)\hat{M}_{\alpha}^{\dagger}\right)$, where $\mathcal{N}$ is a normalization constant. The measurement operators form a proper POVM Nielsen and Chuang (2000), in that

[TABLE]

(i.e. the readouts we have defined here constitute another complete set of measurement outcomes).

It will be useful to take a closer look at the probability density from which the readouts are drawn. Following the procedure we have typically used in the context of optimal paths (OPs) Chantasri et al. (2013); Chantasri and Jordan (2015); Jordan et al. (2015); Lewalle et al. (2017, 2018), we will expand the log of the probability density to $O(dt)$, defining a term

[TABLE]

such that $\wp=e^{C+\mathcal{G}dt+O(dt^{2})}$. We see that up to the two last terms in Eq. (38), the probability density is Gaussian in both readouts, with variances $1/dt$, and means $x\sqrt{\gamma/2}$ and $y\sqrt{\gamma/2}$ for $r_{X}$ and $r_{P}$, respectively. Notice that this corresponds precisely to what we had from the SME, as in (32); the Gaussian form implicit in (38) is in fact key in demonstrating that the form of the SME (23) written in terms of Weiner increments $dW$ is formally suitable for this system.

We also use the Kraus operator to obtain some equations of motion. Consider the exapansion of the Kraus operator itself to $O(dt)$, which reads

[TABLE]

for $|\mathbf{r}|^{2}=r_{X}^{2}+r_{P}^{2}$. We can strip the Gaussian factor $e^{-|\mathbf{r}|^{2}dt/4}$ from the operator for this purpose, since it appears in both the numerator and denominator of the state update expression (11), and thereby cancels off. Consider the following series of approximations, assuming small $dt$:

[TABLE]

which can then be rearranged according to $\rho(t+dt)-\rho(t)\approx dt\>\dot{\rho}$, such that

[TABLE]

This can be expressed in Bloch coordinates by

[TABLE]

It is then straightforward to make the substitutions $r_{X}=x\sqrt{\gamma/2}+\xi_{X}$ and $r_{P}=y\sqrt{\gamma/2}+\xi_{P}$ (32), and see that these equations from the Kraus operator approach are identical to the Stratonovich equations (28) obtained by conversion from the SME approach; this relationship between a Kraus operator based on Bayesian logic, and the SME, is consistent with previous results for this particular measurement Jordan et al. (2015), and other types of continuous qubit measurements leading to diffusive SQTs Gambetta et al. (2008); Chantasri and Jordan (2015); Korotkov (2016); Lewalle et al. (2017).

Simulations can be generated by applying the state update rule (11) with $\hat{M}_{r}\rightarrow\hat{M}_{\alpha}$, with a pair of readouts drawn from Gaussians of means and variances described above, at each timestep. The resulting stochastic trajectories diffuse as expected, and recreate the required decay dynamics on average, as shown in Fig. 5(b).

IV.3 Generalizations

We consider the addition of a Rabi drive to the qubit (i.e. we now discuss additional tones inducing a unitary rotation in the Bloch sphere) by the addition of a Hamiltonian term $i[\rho,\hat{H}]$ to the SME ($\hbar\rightarrow 1$), or a corresponding operator $\hat{U}=e^{-i\hat{H}dt}$ to the measurement scheme with the Kraus operator (where the resulting equations of motion are insensitive to the order of operations, since they are only to $O(dt)$). Without loss of generality, we use $\hat{H}=\delta\hat{\sigma}_{z}/2+\Omega\hat{\sigma}_{y}/2$, where we have denoted the detuning $\delta=\omega_{\text{qb}}-\omega_{\text{dr}}$, with $\omega_{\text{dr}}$ the frequency of the tone555Note that this description is associated with a frame rotating at frequency $\omega_{\text{dr}}$, or equivalently the interaction picture with respect to $\hat{H}_{\text{frame}}=\omega_{\text{dr}}\hat{\sigma}_{z}/2$. In the fixed frame, the qubit Hamiltonian in presence of the drive reads $\hat{H}(t)=\omega_{\text{qb}}\hat{\sigma}_{z}/2+\Omega(i\hat{\sigma}_{-}e^{i\omega_{\text{dr}}t}-i\hat{\sigma}^{+}e^{-i\omega_{\text{dr}}t})$.. Such a tone induces a rotation around an axis tilted by an angle $\text{arctan}(\Omega/\delta)$ with respect to the $z$–axis. Note that the assumption in our derivation has been that only photons emitted by the qubit enter the transmission line which leads to the measurement apparatus; the simplest way to imagine engineering a system such that this remains valid with the Rabi drive on, is that the drive is being implemented by a tone which is off–resonant with the qubit/cavity/transmission line, such that the qubit photons couple to the output leading the measurement device only, and the drive photons couple to their own output only. As discussed earlier, this assumption also requires us to have the cavity resonance far from any of the Mollow triplet peaks, which are centered around $\omega_{\text{dr}}$ and $\omega_{\text{dr}}\pm\Omega_{\text{eff}}$, where $\Omega_{\text{eff}}=\sqrt{\Omega^{2}+\delta^{2}}$ is the generalized Rabi frequency; this regime and assumption is necessary if we want to treat the form of the decay channel as being unaffected by the drive. Drives of the type we have discussed apply generically in “resonance fluoresence” scenarios Kimble and Mandel (1976); Naghiloo et al. (2017); Quijandría et al. (2018), as well as any other situation in which additional tones are present in qubit’s cavity (e.g. to implement additional measurements Ficheux et al. (2018)). The situation we have described here is illustrated in Fig. 1(c).

Note that it is possible to generalize this heterodyne measurement by choosing the phase $\theta$ of the LO. The phase $\theta$ is a relative phase between the signal and LO, so it is equivalent to think of a phase plate having been put in the signal line instead of the LO such that $\hat{a}^{\dagger}\rightarrow e^{-i\theta}\hat{a}^{\dagger}$ in (34), with interference against a fixed pump. Mathematically, we can then assign readouts according to

[TABLE]

which leads to

[TABLE]

The operators for the SME which match the pair of observables we infer from the means of $\mathcal{G}$ are

[TABLE]

We see that changing the phase $\theta$ between the signal and LO effectively rotates the quadrature pair we measure. We have here used notation such that $r_{I}=r_{X}$ and $r_{Q}=r_{P}$ for the choice $\theta=0$. The relationship between the Kraus operator equations of motion and SME equations of motion (Itô or Stratonovich) which we found in the $\theta=0$ case above, hold for arbitrary $\theta$.

We have reviewed the most basic features of an idealized heterodyne measurement. For a more advanced treatment of this system, refer to Jordan et al. (2015); we will now turn our attention to applying the framework we have just developed to homodyne measurement.

V Single–Qubit Homodyne Fluorescence Trajectories

Several experiments Naghiloo et al. (2017, 2016); Tan et al. (2017); Naghiloo et al. (2020) and some theory Bolund and Mølmer (2014) have been published about homodyne fluorescence measurement; we will develop our theory examples here far enough to compare them directly with the simplest experimental results.

V.1 Kraus Operator and Measurement Dynamics

Homodyne detection again involves interfering our signal with a strong LO. Practically, instead of amplifying both quadratures of the resulting signal as in heterodyne detection, homodyne detection involves amplifying one quadrature and de-amplifying the other Tyc and Sanders (2004) (“phase-sensitive” amplification). This procedure amounts to squeezing out the quadrature that isn’t measured, such that in the limit of ideal squeezing we project our signal onto a single quadrature’s eigenstate, instead of onto a coherent state Wiseman (1996); Wiseman and Milburn (2010). This yields a single readout signal, rather than the pair which arise in the heterodyne case. We will follow the same recipe as in the heterodyne case, except that we project onto a final state $\left|X\right\rangle$ (the eigenstate of the $\hat{X}=(\hat{a}^{\dagger}+\hat{a})/\sqrt{2}$ operator in the quadrature space), instead of the coherent state $\left|\alpha\right\rangle$. Again, those unfamiliar with this phase space terminology may wish to consult e.g. Ref. Silberhorn (2007). For dimensionless $X$, recall that we have the following solutions to the quantum harmonic oscillator, which models the field mode:

[TABLE]

Projecting onto the general fluorescence operator, and suppressing the factors $\pi^{-\tfrac{1}{4}}$ on all terms, we get

[TABLE]

Then using a readout substituted in according to

[TABLE]

we find that the POVM is normalized, i.e.

[TABLE]

The relationship between $X$ and $r$ is again set based on comparing the readout statistics with LO’s shot noise, as discussed in appendix A. Those readout statistics can be readily understood from the expression

[TABLE]

which again comes from expanding the logarithm of $\text{tr}(\hat{M}_{x}\rho\hat{M}_{x}^{\dagger})$. We infer that projecting onto $\left|X\right\rangle$ in the photon space leads to a signal related to $x$ in the qubit space, since $r$ has a mean $\sqrt{\gamma}x$, and variance $1/dt$.

As before, we may take $\hat{a}^{\dagger}\rightarrow\hat{a}^{\dagger}e^{-i\theta}$ to generalize the choice of measured quadrature, yielding an operator

[TABLE]

which still generates a proper POVM. Expanding the log–probability density for the readout gives us

[TABLE]

thus the mean of the Gaussian in $r$ matches the signal given by $\hat{L}+\hat{L}^{\dagger}=\sqrt{\gamma}(\hat{\sigma}_{x}\cos\theta-\hat{\sigma}_{y}\sin\theta)$ for the SME operator $\hat{L}=\sqrt{\gamma}e^{-i\theta}\hat{\sigma}_{-}$.

We proceed to find the equations of motion. Note that we can approximate $\hat{M}_{x}$ as we did $\hat{M}_{\alpha}$ (39), such that

[TABLE]

Then by the logic of (40) and (41) we may derive equations of motion in terms of Bloch coordinates

[TABLE]

We have again used a Rabi drive characterized by $\hat{H}=\Omega\hat{\sigma}_{y}/2+\delta\hat{\sigma}_{z}/2$, or $\hat{U}\approx\openone-i\Omega\hat{\sigma}_{y}dt/2-i\delta\hat{\sigma}_{z}dt/2$. As above, these equations are consistent with those derived from the SME (23), using $\hat{L}=\sqrt{\gamma}e^{-i\theta}\hat{\sigma}_{-}$, provided the SME output is correctly interpreted as an Itô equation, whose Stratonovich form then matches the above exactly. Simulated trajectories for the case $\theta=0$ and $\Omega=0$ are shown in Fig. 5(c), and demonstrate good agreement with expectations, as in the previous cases.

V.2 Inefficient Measurements

Inefficient measurement is easily included in the SME (23), and is completely described by the dimensionless parameter $\eta\in[0,1]$. In the Kraus operator picture, we must modify the amplitude of the signal going into the measurement apparatus; we will find a case intermediate between perfect measurements (11) and no measurement (20), reflecting that some fraction of the information is lost rather than collected. A straightforward way to represent this is with an unbalanced beamsplitter placed in front of our (still otherwise ideal) measurement device, as shown in Fig. 1(d). If $\hat{a}^{\dagger}$ creates a photon in the emitted field mode, the beamsplitter transforms it according to

[TABLE]

where the surviving signal $\hat{a}_{s}^{\dagger}$ goes to the detector with probability $\eta$, but the information in channel $\hat{a}_{\ell}^{\dagger}$ is lost with probability $1-\eta$, and outcomes (all of which could have occurred) in the latter channel must be traced out. We will do the trace of the lost channel in the Fock basis for simplicity (a sum of two terms is simpler than an integral over a continuous homodyne or heterodyne readout, although averaging over any complete set of hypothetical measurement outcomes is technically correct). The scheme we are describing, for homodyne detection with efficiency $\eta$, can be implemented with a pair of operators

[TABLE]

for Fock states $j=0,1$ in the lost mode, i.e.

[TABLE]

with a state update rule

[TABLE]

The measured homodyne signal is computed according to projection onto the states $\left|X\right\rangle$ exactly as above, and a drive could be added with unitaries in the same manner as above. The new operators (56) again denote a well–defined measurement, in that they form POVM elements, i.e.

[TABLE]

We find the same agreement between the expansion of the state update (58) to $O(dt)$, and the SME with finite $\eta$ (converted to its Stratonovich form), as in every case discussed. Thus the description of $\eta$ supposed by Fig. 1(d) and (55) is entirely equivalent to the description implicit in the SME, and clarifies the meaning of measurement “inefficiency”.

This picture of inefficiency is also readily connected to scenarios in which several observers simultaneously make measurements, and each gets only partial information Jacobs and Steck (2006); Dziarmaga et al. (2004); Harrington et al. (2019). One can imagine that an observer lives at each output of the beamsplitter in Fig. 1(d), each recieving some proportion of the information about the qubit carried by the decay process as they make measurements. If they do not share their results, each will have a different estimate of the qubit’s evolution conditioned on their partial information, and tracing out over the other observer’s measurement record which they do not have access to. Either of their estimates could be compared to some hypothetical “true” evolution which an observer able to access all the relevant measurement records could compute. In practice, it is effectively impossible to have a perfectly efficient measurement in any experiment, and some information is always irretrievably lost to the environment through any channel from which the primary system is not perfectly isolated (generically, this is “decoherence”). The methods we have presented here can readily be adapted to the kind of multiple–observer situation we have just described; this includes situations which involve both jumps and diffusion, due to different observers making different types of measurements (see e.g. appendix B of Lewalle et al. (2019a), or Kuramochi et al. (2013)). Such scenarios have recently been fruitfully investigated in the context of quantum state smoothing Guevara and Wiseman (2015); Chantasri et al. (2019); Guevara and Wiseman 666Quantum state smoothing is closely related to quantum trajectories; SQTs, as we have presented them in the present text, are a form of “quantum filtering” which goes forward in time; in other words, we here only use the measurement record from the system’s past to estimate a qubit’s state. In the event of an inefficient measurement, quantum state smoothing often allows for a more pure estimate of the system’s state to be made at some time, by using the measurement record both before and after the time at which the state is estimated..

We perform simulations which include measurement inefficiency, which are shown in Fig. 7, and discussed further in connection with the “optimal path” techniques we develop shortly. Measurement inefficiency leads to decay that is qualitatively the same as in the ideal case discussed in Fig. 5, except that instead of trajectories being restricted to pure states on the surface of the Bloch sphere, as in the $\eta=1$ case, they instead move stochastically on the surface of an ellipsoid which contracts towards $\left|g\right\rangle$ over time as information is lost in the $\eta<1$ case. Qualitative agreement between these simulated results shown in Fig. 7 and those obtained in experiment for either the homodyne Naghiloo et al. (2016) or heterodyne Campagne-Ibarcq et al. (2016a) detection can be verified at a glance, and a quantitative understanding of this will be developed shortly.

As we have now successfully adapted and extended our presentation of basic methods for pertaining to quantum trajectories to the homodyne detection case, and established that they behave correctly, we can proceed by extending our analysis of this system into new examples which can introduce and highlight particular topics in the recent literature.

VI Special Topics and Further Examples

We will focus on connections to two areas, using homodyne fluorescence detection as our example of choice; first we describe how this example relates to recent work about the arrow of time in quantum measurement, which connects to work on fluctuation theorems for quantum trajectories, and the growing area of quantum thermodynamics more generally; second we will describe how “most–likely paths” can be derived from diffusive quantum trajectory dynamics using a variational principle.

VI.1 Time reversal symmetry and the arrow of time

How does an arrow of time emerge from microscopically time–reversible physical laws? The issue has been raised in the context of continuous quantum measurements Dressel et al. (2017); Manikandan and Jordan (2019); Manikandan et al. (2019); Harrington et al. (2019), and applies more broadly across many disciplines within physics Hawking (1985); Maccone (2009); Lebowitz (1993). In the quantum measurement case, one could pose this question as a game; a quantum trajectory is shown like a movie, forward and backward, and the goal of the game is to infer the direction in which the movie was originally recorded. We will find that the equations of motion are time–symmetric, (e.g. as in Hamiltonian dynamics), such that both the forward and backward movies both depict legitimate dynamics; this is not the whole story however, as the backward evolution (i.e. “wavefunction uncollapse”) does not necessarily occur with the same probability Jordan and Korotkov (2010); Korotkov and Jordan (2006). This leads to a natural discriminator for the arrow of time in terms of the probabilities of occurrence of forward and backward trajectories of the monitored quantum system, as developed in Refs. Dressel et al. (2017); Manikandan and Jordan (2019). Assuming no prior bias, we could use the measurement record as an additional tool to improve our inference about the direction in which the quantum state dynamics is originally recorded Dressel et al. (2017) (by analogy, the sound track for a movie could help us understand in which direction it is meant to run). Such an approach is fundamentally connected to the time–symmetry of underlying dynamical equations describing the measurement, and connects to the arrow of time analysis pertinent to the thermodynamics of small systems Chernyak et al. (2006); Manikandan et al. (2019).

Fluorescence appears to exhibit a clear arrow of time, and therefore the time–reversibility of continuously monitored fluorescence dynamics may seem rather surprising. For this reason, here we detail the time symmetry analysis of dynamical equations (54) which describe homodyne measurement of fluorescence, using the approach presented in Ref. Manikandan and Jordan (2019), wherein similar and detailed analysis was performed for the heterodyne case. The time–reversed dynamics can be considered as a legitimate measurement dynamics, starting from the time–reversed final state $\Theta|\psi_{f}\rangle$, evolving through the time–reversed counterpart of the forward sequence of states, back to the time–reversed initial state $\Theta|\psi_{i}\rangle$. The measurement operators of the backward dynamics are related to the forward dynamics by a Hermitian conjugate operation, i.e. $\hat{M}_{B}=\hat{M}_{F}^{\dagger}$; therefore the dynamical equations which describe the backward dynamics are also similar to the retrodicted dynamical equations Tan et al. (2015), but starting from the time–reversed final state. We may write the retrodicted dynamical equations corresponding to a homodyne measurement, where the quantum state is updated by

[TABLE]

We have parameterized the single–qubit density matrix \underaccent{\bar}{\rho} with Bloch coordinates according to

[TABLE]

Using the form of measurement operators given in Eq. (51) The dynamical equations now take the form,

[TABLE]

Note that the retrodicted equations under the time-reversal operation, \underaccent{\bar}{x}\rightarrow-x,~{}\underaccent{\bar}{y}\rightarrow-y,~{}\underaccent{\bar}{z}\rightarrow-z, and \underaccent{\bar}{t}\rightarrow T-t (i.e., \underaccent{\bar}{dt}\rightarrow-dt) looks exactly like the forward dynamical equations (54), demonstrating their time–reversal invariance; we have eliminated the drive characterized by $\Omega$ and $\delta$ in the equations above for brevity, but including it does not affect the result. Time reversal symmetry of the dynamical equations suggests that the forward dynamics and the reverse dynamics both represent a physical quantum trajectory on the Bloch sphere. Given the measurement record, one can associate a probability each to the forward and backward trajectories, which can be used to infer an arrow of time for the measurement dynamics, and subsequently characterize the irreversibility of homodyne measurement of fluorescence using the associated fluctuation theorems Dressel et al. (2017); Manikandan and Jordan (2019); Manikandan et al. (2019).

We can expand on this story somewhat by noting that our diffusive trajectories in Fig. 5(b,c) do not diffuse monotonically downward from $\left|e\right\rangle$ towards $\left|g\right\rangle$. This suggests that the measurement process can actually cause the probability of the qubit being found in the more energetic of its states to rise in some realizations. While the average decays monotonically, fluctuations make this question of an arrow of time non–trivial in individual realizations; the probability of sustained re-excitation over a long period is low, but estimates of the arrow of time using only a short window of the evolution cannot necessarily be made with high confidence; rare events can decieve, and something resembling a “wavefunction uncollapse” is not merely hypothetical in this system. Such behavior has been noted in the literature Bolund and Mølmer (2014); while it may be initially intuitively challenging, this effect is perfectly correct, and reflects the nature of the information we get about the field when we make a weak ($dt\ll T_{1}$) quadrature measurement, and its backaction. A truly detailed description of the thermodynamics of quantum measurements or trajectories falls beyond our present scope, but is a fascinating area related to the questions we have discussed here, and enjoying increased recent research interest Cottet et al. (2017); Elouard et al. (2017a, b); Naghiloo et al. (2020); Elouard and Jordan (2018); Masuyama et al. (2018); Monsel et al. (2018, 2019); Mohammady et al. (2019). We encourage the curious reader to explore further.

VI.2 A Variational Principle for Quantum Trajectories

Optimal paths (OPs) Chantasri et al. (2013); Chantasri and Jordan (2015); Jordan et al. (2015); Areeya Chantasri (2016); Lewalle et al. (2017, 2018) have recently been used to elucidate a variety of quantum trajectory phenomena. We will here give a brief overview of their derivation, using the CDJ (Chantasri/Dressel/Jordan) path integral Chantasri et al. (2013), and then apply the formalism to the homodyne fluorescence examples we have developed above.

VI.2.1 Derivation of Optimal Paths

OPs can be understood as the path extremizing the probability to get from one given quantum state $\mathbf{q}_{i}$ to another $\mathbf{q}_{f}$ in a particular time interval, under the dynamics due to backaction from the continuous weak quantum measurement. The vector $\mathbf{q}$ parameterizes the quantum state, and here denotes coordinates on the Bloch sphere. Typically OPs will be most–likely paths (MLPs), which maximize the probability of the measurement record connecting the given boundary conditions according to an action–extremization principle. OPs should be confused neither with a globally most–likely path (i.e. the particular MLP post–selected on the most likely final state after a given time interval), or with an average path. Details about numerical procedures to extract an MLP from data, which corresponds to the theory we are about to develop, can be found in appendix B.

To begin, we explain how the path probability can be written in terms of an effective action, which can then be extremized according to a variational principle. We may write an expression for the probability of a quantum trajectory, which moves from $\mathbf{q}_{i}$ to $\mathbf{q}_{f}$ through a discrete sequence of measurements as

[TABLE]

The $\delta$–functions at the initial and final points impose the boundary conditions. The indices $k$ run over time, such that if $\rho_{k}=\rho(t)$, then $\rho(t+dt)=\rho_{k+1}$ and so on. The stochastic element of the dynamics arises in drawing the readout from the probability density determined by the denominator of the state update expression (11), i.e. $\wp(\mathbf{r}_{k}|\rho_{k})\propto\text{tr}(\hat{M}_{\mathbf{r}_{k}}\rho_{k}\hat{M}_{\mathbf{r}_{k}}^{\dagger})$, where $\hat{M}_{\mathbf{r}_{k}}$ could generically be any Kraus operator describing a weak measurement. We describe the deterministic update of the quantum state given the stochastic readout $\mathbf{r}_{k}$ according to $\wp(\mathbf{q}_{k+1}|\mathbf{q}_{k},\mathbf{r}_{k})=\delta(\mathbf{q}_{k+1}-\mathbf{q}_{k}-dt\>\mathcal{F}[\mathbf{q}_{k},\mathbf{r}_{k}])$, where $\dot{\mathbf{q}}=\mathcal{F}[\mathbf{q},\mathbf{r}]$ is an equation of motion, e.g. like (54). Recall that a $\delta$–function may be written $\delta(\mathbf{q})=(2\pi i)^{-\text{dim}(\mathbf{q})}\int_{-i\infty}^{i\infty}d\mathbf{p}\>\exp\left[-\mathbf{p}\cdot\mathbf{q}\right]$, where $\text{dim}(\mathbf{q})$ is the dimension of $\mathbf{q}$, and $d\mathbf{p}=dp_{1}\;dp_{2}\;...\;dp_{\text{dim}(\mathbf{q})}$. We apply this identity to all $\delta$–functions in (63), such that in the time–continuum limit we have a Feynman–like path integral, in which we effectively sum over all possible quantum trajectories which obey the given boundary conditions, i.e.

[TABLE]

for $\mathcal{N}=(2\pi i)^{-(n+2)\cdot\text{dim}(\mathbf{q})}$, and where $\mathcal{D}[\mathbf{p}]$ arises from the infinite product of the $d\mathbf{p}_{k}(2\pi i)^{-\text{dim}(\mathbf{q}_{k})}$. We use the shorthand $\mathcal{B}=-\mathbf{p}_{-1}\cdot(\mathbf{q}_{0}-\mathbf{q}_{i})-\mathbf{p}_{n}\cdot(\mathbf{q}_{n}-\mathbf{q}_{f})$ for the boundary terms, and the shorthand $\mathcal{G}$ for the expansion to $O(dt)$ of the log–probability for the readouts $\ln\wp(\mathbf{r}|\mathbf{q})$ (see e.g. (52)). This relates a trajectory probability $\mathcal{P}$ to a “stochastic action” $\mathcal{S}$. That action is expressed in terms of a Hamiltonian $\mathcal{H}=\mathbf{p}\cdot\mathcal{F}+\mathcal{G}$.

We can then see that extremizing the probability corresponds to extremizing $\mathcal{S}$. Action extremization in this case can be expressed much the same way as it is in classical mechanics, such that

[TABLE]

which indicates that the OPs obey the equations

[TABLE]

These are Hamilton’s usual equations in the Bloch coordinates $\mathbf{q}$ and conjugate variables $\mathbf{p}$, plus an additional equation which stipulates that the stochastic readouts be optimized, leading to some smooth $\mathbf{r}^{\star}$ instead of the stochastic $\mathbf{r}$. The OPs are thus smooth curves, and are solutions to a Hamiltonian dynamical system of ordinary, rather than stochastic, differential equations. The OPs are themselves possible quantum trajectories (since $\dot{\mathbf{q}}=\mathcal{F}$ is preserved, by construction), even though the optimal readouts $\mathbf{r}^{\star}$ are smooth functions of time, rather than the stochastic readouts which occur in individual runs of an experiment. The Hamiltonian structure of the OP dynamics implies that, absent any explicitly time–dependent parameters, there is a conserved “stochastic energy” $E=\mathcal{H}$ associated with an OP.

The “momenta” $\mathbf{p}$ conjugate to the generalized Bloch coordinates $\mathbf{q}$, which arise in this optimization process, warrant further attention. The $\mathbf{p}$ are not directly measurable, but play a substantive mathematical role, in that they effectively generate displacements in the quantum state $\mathbf{q}$ according to the optimal measurement record $\mathbf{r}^{\star}$. They appear as the variables conjugate to $\mathbf{q}$ in the Fourier representation of the $\delta$–functions, and could be understood and Lagrange multipliers in the optimization process. They are perhaps best understood in analogy with classical optics, however: just as we have expressed a diffusive process in terms of an action, the underlying wave process in classical optics can be expressed as an action; extremization in the latter scenario leads to a ray description. Our OPs are related to the underlying diffusive process described above in much the same way that a ray description of light is related to the underlying wave optics process777We are grateful for comments by Prof. Miguel Alonso which helped us to clarify this point in our own thinking.. Diffusive SQTs arrive at final states with a variety of probabilities (and the action representing this is real), whereas optical paths leading to different positions arrive with different phases, exhibiting interference (and the subsequent action in this case is imaginary). Continuing with this analogy, we could imagine that rays leave a source in different directions, i.e. with different wave–vectors (OPs leave their initial state $\mathbf{q}_{i}$ with a range of initial momenta $\mathbf{p}_{i}$); we may pick out a particular ray or subset of rays by choosing a particular initial wave–vector, or by choosing a final position to which they connect (a particular OP can be selected by choosing a value of $\mathbf{p}_{i}$, or by choosing a $\mathbf{q}_{f}$ at a later time). Thus, the degree of freedom in choosing an initial momentum $\mathbf{p}_{i}$ is the same degree of freedom which allows for post–selections to many $\mathbf{q}_{f}$; the mapping between the two is not necessarily one–to–one (see our work on “multipaths” for clarification of this point Lewalle et al. (2017); Naghiloo et al. (2017); Lewalle et al. (2018)). In a phase space with $N$ coordinates $\mathbf{q}$ and $N$ momenta $\mathbf{p}$ (a $2N$–dimensional phase space), the set of paths evolving from all $\mathbf{p}_{i}$ with a fixed $\mathbf{q}_{i}$ defines an $N$–dimensional Lagrangian manifold (LM) within the phase space; such a manifold can be understood as containing all of the possible dynamics originating from the state $\mathbf{q}_{i}$ in the OP description, i.e. such an LM explores the full space of optimal readouts $\mathbf{r}^{\star}$, just as the underlying diffusive process may explore the full space of stochastic measurement records888The LM in question has primarily been used in the context of multipath dynamics Lewalle et al. (2017); Naghiloo et al. (2017); Lewalle et al. (2018); there the main concern is whether the projection of the time–evolved LM out of the full OP phase space down into the $\mathbf{q}$–space of final quantum states is one–to–one (a single MLP connects the initial state to the chosen final state) or many–to–one (in which case many OPs may connect the boundary conditions, typically corresponding to different clusters of SQTs in the post–selected distribution).. Our focus now will be on applying this formalism to our homodyne fluorescence measurement; some work in this vein, albeit with a different emphasis, appears in Ref. Naghiloo et al. (2017).

VI.2.2 Optimal Paths for Homodyne Fluorescence Trajectories

Notice that the system of equations (54) can be simplified straightforwardly; $\dot{y}=0$ if $y=0$ and $\theta=0$, and then all the dynamics are in the $xz$–plane of the Bloch sphere. It is then easy to write down a stochastic Hamiltonian

[TABLE]

for the OP dynamics using formulas we have already derived above (with $\Omega=0=\delta$, and $\eta=1$) to describe the ideal measurement, its backaction, and statistics999Note that we have derived both $\mathcal{F}$ and $\mathcal{G}$ using regular calculus, and are thereby effectively using the Stratonovich form of $\mathcal{F}$, not the Itô form which arises directly from the SME (23). Using a form of $\mathcal{F}$ which does not transform according to regular calculus would prevent us from performing our OP analysis using typical approaches of classical mechanics (e.g. canonical transformations), which we find quite undesireable.. The optimal readout obeys $\partial_{r}\mathcal{H}_{hom}^{xz}|_{r^{\star}}=0$, which we solve to obtain

[TABLE]

we see that we have the signal $\sqrt{\gamma}x$, plus some additional terms which depend on the conjugate momenta $p_{x}$ and $p_{z}$, which implement the optimized effect of the noise (as discussed above).

We can simplify the equations even more. Consider a change to polar coordinates according to the canonical transformation

[TABLE]

which preserves the Poisson brackets between all of the pairs of conjugate variables. Then we see that for the choice $R=1$ and $p_{R}=0$, we have $\dot{R}=\partial_{p_{R}}\mathcal{H}_{hom}^{\star,R\vartheta}=0$, meaning that we can look at dynamics purely on the great circle of the Bloch sphere at $R=1$ (pure states), where states are parameterized entirely by a single coordinate $\vartheta$ (with $\vartheta=0\leftrightarrow\left|e\right\rangle$ and $\vartheta=\pi\leftrightarrow\left|g\right\rangle$). Making this transformation, and subtituting in the optimal readout such that $\mathcal{H}^{\star}=\mathcal{H}|_{r=r^{\star}}$, we obtain the Hamiltonian

[TABLE]

which generates the OPs in the simplest case we can consider for this system.

The phase space for this Hamiltonian is plotted in Fig. 6, along with the time–derivative of the stochastic action $\mathcal{S}$ which is extremized by the OP dynamics (effectively, $\dot{\mathcal{S}}$ gives an approximate representation of the probability cost involved with traversing certain regions of the OP phase space). A careful reading of these plots can provide an insightful overview of the system dynamics. First, we can immediately infer a rule of thumb: OPs with higher stochastic energy generically correspond to events which occur with lower probabilities (this is true to the extent that regions of large $E$ correspond to regions of more–negative $\dot{\mathcal{S}}$). Secondly, we see that all paths in the OP phase space eventually approach $\left|g\right\rangle$ ($\vartheta=\pi$) in the long–time limit, as we expect they must; there are possibilities for this to occur in either direction around the Bloch sphere with some probabilities, but these pure–state OPs never cross through $\vartheta=\pi$. The uni–directionality of the flow towards $\left|g\right\rangle$ after $t\gg T_{1}$ reflects our intuition that there should be a statistical arrow of time in the measurement–induced dynamics, as is discussed above and in detail elsewhere Manikandan and Jordan (2019); Manikandan et al. (2019). A particular point in the phase space is worthy of further attention; the unstable fixed point at $\vartheta=0$ and $p=0$ describes an OP which is stationary at $\left|e\right\rangle$ for all time; this does not violate our intuition however, since the probability cost involved with sitting at that point is greater than for sitting very close to $\left|g\right\rangle$, such that it is still virtually impossible to post–select on a state still at $\left|e\right\rangle$ after $t\gg T_{1}$. As noted before, it is possible for paths to start near the ground state, and re–excite, passing through $\left|e\right\rangle$ before asymptotically approaching $\left|g\right\rangle$ again around the other side of the Bloch sphere; while such behavior corresponds to relatively rare events (a relatively low–probability post–selection is required), the possibility of such events is readily visible in the OP phase portrait101010For example, starting at $\vartheta_{0}=\pi-1$, a globally more–likely path which we would regard as typical might arise from post–selecting on $\vartheta_{f}=\pi-0.01$ ($\approx$ ground) at some later time, picking out a solution from X+. However, it is possible to post–select on some state like $\vartheta_{f}=-\pi+0.01$ (also $\approx$ ground) instead, selecting a path from X- or Y-; this reveals the possibility of a much rarer set of events, corresponding to an OP which circles back through the excited state before decaying towards ground, from the opposite direction as compared with the more typical set of paths. Post–selections drawing out these dynamics can generically select OPs from regions X or Y of the phase portrait, as detailed in Fig. 6; those in region Z can only partially re-excite, before turning around to decay.. Further details about OPs this system, for the case of $\Omega\neq 0$, can be found in Naghiloo et al. (2017), and more detailed investigations of the corresponding heterodyne cases can be found in Jordan et al. (2015).

VI.3 Optimal Paths for Inefficient Measurements, and Connections to Experiments

We develop a final example; in order to compare OPs directly to experimental results, we need to introduce measurement inefficiency. In the homodyne case with $\theta=0$ and $\Omega=0=\delta$, corresponding to the SME readout $r=\sqrt{\eta\gamma}\>x+\xi$, and the equations of motion from either expanding (58) to $O(dt)$, or converting the requisite SME from Itô to Stratonovich, are

[TABLE]

We immediately see that for this choice of measured quadrature, the $y$ component of the dynamics can be eliminated with the choice $y=0$, leaving only dynamics in the $xz$–plane of the Bloch sphere; we will assume $y=0$ for the remainder of this section. These assumptions, with imperfect $\eta$, give us the simplest version of this system that can be compared directly with existing experiments. The last piece we need is an understanding of the probability density function from which the readouts are drawn; using the same methods as above, we find

[TABLE]

Thus we see that simulations involve repeated state updates as per (58), with readouts drawn at each step from a Gaussian of mean $x\sqrt{\eta\gamma}$ and with variance $dt^{-1}$. Optimal paths are derived from a stochastic Hamiltonian

[TABLE]

where $f_{x}$ and $f_{z}$ are the RHS of (71a) and (71c), respectively. Our aim below will be to elucidate the basic aspects system dynamics, and show that our simulations and OPs match relevant results in the experimental literature.

A particularly important feature of the dynamics under homodyne fluorescence detection (absent a Rabi drive or other dynamics) is that all trajectories are constrained to an ellipse in the Bloch sphere at any given time Naghiloo et al. (2016); Tan et al. (2017) (and a similar ellipsoid is apparent in the heterodyne case Campagne-Ibarcq et al. (2016b)). The functional form of these ellipses has been derived in the literature Tan et al. (2017), and follows

[TABLE]

for the time–dependent function

[TABLE]

where $u_{0}$ is set by the initial state according to

[TABLE]

For example, with the initial state $\left|e\right\rangle$, we have $u_{0}=1$ at $t=0$, and at any time $t>0$ all possible trajectories evolving from $\left|e\right\rangle$ under dynamics from the inefficient homodyne measurement can be found on the ellipse (74) (as a function of $x$); this ellipse is initially the great circle bounding the $xz$–plane of the Bloch sphere, and decays towards the ground state according to the time dependence (75). We develop this example in Fig. 7. In the left four panels (a–d) we show the density of simulated trajectories originating at $\left|e\right\rangle$ after different evolution times; we find essentially perfect agreement between the histograms of these simulated trajectory densities, and the analytic curves (74) known from the literature. We stress that in departure from many other quantum measurement scenarios, there are final states on which it is impossible to post–select in the present system; typically some states appear very rarely in the dynamics, but here large regions of the Bloch sphere are forbidden entirely.

We conclude by demonstrating the connection between the Lagrangian manifold from the OP phase space we described above, and the ellipses we have just described. The relative simplicity of the dynamics under inefficient homodyne fluorescence measurement make this an ideal example with which to illustrate the concepts discussed above. In Fig. 7(e) we show the projection of the LM originating at $\left|e\right\rangle$ into the $xz$–plane of the Bloch sphere (i.e. we evolve the OP equations sampled across the initial LM, and then flatten the two–dimensional manifold, which lives in the four–dimensional phase space, into a plot that appears in the coordinates $x$ and $z$ only, at selected times). We then see that we have exact agreement between the LM and the analytic curves (74), consistent with the fact that the OPs are themselves possible quantum trajectories. This reinforces our statements about the consistency between the methods reviewed here, and the broader literature on continuous monitoring of fluorescence, but also serves to illustrate the role of the initial momenta $\mathbf{p}_{i}$ in the OP formalism. Choosing a particular $\mathbf{p}_{i}$ selects particular boundary conditions from the possible multitude, and the complete set of $\mathbf{p}_{i}$ contained in the LM index a complete set of possibilities for the OPs originating at a particular state. While the LM in question may at first seem a somewhat abstract mathematical object, we are here able to highlight its physical character.

VII Closing Remarks

We have given an overview of many useful methods and insights that arise from considering continuous quantum measurement, emphasizing examples in which we track a quantum emitter’s state by gathering and measuring spontaneously–emitted photons. We have focused on a Kraus operator approach to this problem, most similar to that developed in Ref. Jordan et al. (2015), and made connections to a corresponding stochastic master equation description throughout. Many of the issues which arise in treating this particular type of system are common to stochastic quantum trajectories in general, and we have consequently addressed many of the important principles and typical problems one needs to become aware of when entering this research area. We have also been able to use the fluorescence examples above to offer accessible illustrations and introductions to selected advanced topics of contemporary interest; for example, we have been able to make comments to help the interested reader engage with work on the arrow of time in quantum trajectories, or understand how to generate and interpret trajectories which follow an optimal measurement record.

We can take a larger view of the processes we have described. We are accustomed to talking about the fluorescence process in terms of the emission of individual photons at particular times, from a sudden jump in the qubit state, i.e. we typically discuss fluorescence in language which lends itself naturally to the photodetection case shown in Fig. 5(a). This notion of a jump and photon emerging at a particular time were mathematically enforced in the framework presented above by choosing outcomes in the Fock basis of the field, and with perfect measurement efficiency, we can say that we have collected complete information about the output mode. The other measurements we have discussed, are however, just as “complete” (in the sense that we have a POVM, and can ascribe a pure state to the qubit at any time in the $\eta=1$ case; even in the $\eta<1$ case, we are able to assign a state to the qubit which is consistent with a subsequent tomographic verification Campagne-Ibarcq et al. (2016b); Naghiloo et al. (2016)). The trajectories generated by the heterodyne and homodyne measurement schemes do not readily admit interpretation in terms of a photon emerging at any particular time; there is no single point along a trajectory in Fig. 5(b or c) to which we can point and say “this is when the photon was emitted”. However, we can still ascribe a stochastic evolution to the qubit state, which agrees with any sensible check we know how to perform, and that evolution evidently depends on how or what type of environmental information was collected; the dynamics on average reflect the same decay statistics we are used to regardless of the measurements (which is always the case in such continuous quantum measurement problems, and more or less the only assumption we made to formulate the model we have used). In other words, we the observer can dynamically assign our qubit a state111111One could say that we, the observer, are updating our probabalistic prediction about outcomes of a future measurement; this “best guess” about possible future outcomes is made based on some preparation procedure, past measurement record, and the rules we typically refer to as “quantum mechanics”. which may e.g. begin at $\left|e\right\rangle$ and eventually wind up in $\left|g\right\rangle$ according to the typical statistics of spontaneous emission; none of the other details of our photodetection story carry over in a simple way, however, to the case of generalized measurements. All of this is related to the fact that the emission process entangles the qubit with the field mode (see e.g. the transition from (12) to (13)); this correlation allows us, as observers, to infer an emitter state based on the field, but also means that the type of information we infer about the emitter depends on what kind of state we project the field into (we must choose what kind of question to ask the field, and this will affect the kind of answer we can subsequently expect). Such issues get to the heart of quantum mechanics, and quantum trajectories generally serve as an excellent point of departure for such discussions about the foundations and interpretations of the theory, and the role of an observer probing a system which is otherwise sufficiently isolated so as to behave “quantumly”.

There are, of course, many related topics which we have not been able to cover at all, but which we hope may become substantially easier for a new reader to digest with the foundation we have developed above. We highlight three in particular, in addition to those of a thermodynamic character Elouard et al. (2017a, b); Naghiloo et al. (2020); Masuyama et al. (2018), pertaining to quantum state smoothing Guevara and Wiseman (2015); Gough (2020); Chantasri et al. (2019); Guevara and Wiseman , or the optimal path methods Chantasri et al. (2013); Chantasri and Jordan (2015); Jordan et al. (2015); Areeya Chantasri (2016); Lewalle et al. (2017, 2018) we have already mentioned. First, some of the main interest in quantum trajectories is geared towards its applications to quantum control; we point out some literature which adapts the types of measurements we have described above to this purpose Wang and Wiseman (2001); Campagne-Ibarcq et al. (2016b). Second, dispersive measurements Murch et al. (2015); Gambetta et al. (2008); Murch et al. (2013); Weber et al. (2014), allowing for measurements of e.g. $\hat{L}\propto\sigma_{z}$, are extremely common in the literature; unlike the measurements we have described here, a qubit in a cavity (or coupled to a resonator) is directly probed with a pulse which reveals information about its state; virtually everything in terms of general approach we have developed above can be carried over to this case however. Such measurements have also been utilized in contemporary feedback control problems Patti et al. (2017); Hacohen-Gourgy et al. (2018); Minev et al. (2019), or in simultaneously weakly measuring multiple non–commuting observables Hacohen-Gourgy et al. (2016); Lewalle et al. (2017, 2018); Atalaya et al. (2018). Combinations of dispersive and fluorescence measurements have been realized experimentally Ficheux et al. (2018). An introduction to a simple model of dispersive measurements, similar in spirit to that above, can be found in the appendix of Ref. Lewalle et al. (2017) and references therein. Third, direct extensions of the measurements we have discussed here to the two–qubit case can serve as a springboard to study measurement–induced entanglement generation between a pair of emitters Mascarenhas et al. (2011); Santos et al. (2012); Lewalle et al. (2019b, a), and the decay of entangled states open to the influence of decay channels and/or measurements Yu and Eberly (2004); Mintert et al. (2005); Viviescas et al. (2010); Mascarenhas et al. (2010); Carvalho and Santos (2011). “Bell state measurements” are essential in recent tests of quantum mechanics and local realism Wiseman (2015). As such, we hope that our present work may help new readers to better digest a wide variety of literature concerning quantum measurements, open quantum systems, and beyond.

Acknowledgements.

We acknowledge funding from NSF grant no. DMR-1809343, and US Army Research Office grant no. W911NF-18-10178. PL acknowledges additional support from the US Department of Education grant No. GR506598 as a GAANN fellow, and thanks the Quantum Information Machines school at École de Physique des Houches for their hospitality during part of this manuscript’s preparation. We are grateful to Joe Murphree for helpful comments. Analyses of the type described in appendix B and Fig. 8 have benefited from numerous conversations with Areeya Chantasri, and some versions of the underlying numerical methods have benefited from work by John Steinmetz and Kurt C. Cylke. Numerical methods underlying Figs. 3, 5, 7, and 8 have been implemented in Python 2.7.

APPENDICES

Appendix A LO Power and Photocurrent Fluctuations

We consider the logic behind the scaling of the readouts in time in the dyne measurements, i.e. we justify the expressions (35) and (48). The argument we offer follows directly from Ref. Jordan et al. (2015). Dyne measurements involve interfering the signal beam (which here contains only zero or one photon) with a strong coherent state LO. Recall that a coherent state

[TABLE]

has a mean photon number

[TABLE]

with fluctuations

[TABLE]

This is a direct result of the fact that a coherent state generates the Poisson statistics for photon arrival times; this is typical for arrival times pertaining to any random process characterized by a constant average rate. A constant underlying rate is, of course, consistent with an assumption that the LO has constant power (on average, up to the quantum fluctuations). If $h\nu$ is the energy per photon at frequency $\nu$ in the LO beam, then the LO power corresponds to the average photon number arriving in a time interval $dt$ according to $\left\langle P\right\rangle=\left\langle N\right\rangle h\nu/dt$, and we have fluctuations $\left\langle N\right\rangle\pm\sqrt{\left\langle N\right\rangle}$ or

[TABLE]

We can talk about the fluctuations either in photon number, or LO power, which are

[TABLE]

Thus we stress that for fixed $\nu$ and $\left\langle P\right\rangle$, the fluctuations in the LO photon number go like $\sqrt{dt}$, and the fluctuations in the LO power go like $1/\sqrt{dt}$. In other words, for fixed average LO power, there are fluctuations in the subsequent photocurrent with variance $1/dt$, and it is precisely these physics which motivate the assignments (35) and (48) in modeling the heterodyne and homodyne measurements, respectively. Those fluctuations are well approximated as Gaussian for large $\left\langle N\right\rangle$ (the Poisson distribution converges to a Gaussian for large numbers). That our formalism leads to agreement with standard tools like the SME, and with experimental data (e.g. as in Fig. 7), further justifies the use of the intuition we provide here, connecting our measurement noise to the shot noise of the measurement LO / amplification pump.

Appendix B Connection Between Optimal Paths and Simulation / Data

Our expectation in deriving most–likely paths (MLPs) is that, in typical cases, they should correspond approximately to a highest–probability peak in the post–selected sub–ensemble of trajectories connecting an initial and final state (see comments below about some of the subtleties, however). Our aim here is to describe the procedure by which we approximate this concept in extracting a MLP from simulated SQTs or data, and show by example that such results are in good correspondence with those given by the CDJ optimization (see Sec. VI.2.1). This has been performed in a variety of cases elsewhere Weber et al. (2014); Naghiloo et al. (2017); Lewalle et al. (2018), and we here describe and perform the requisite analysis for solutions of (70) and the corresponding simulation.

We generically suppose that we are given a simulated ensemble of SQTs $\{{\rho(t)}\}$, initialized at a particular $\rho_{0}$; such trajectories are necessarily sampled over some small but discreet timestep $dt$, consistent with the simulation procedures described in the main body of the text above. We describe the numerical manipulations we perform on such a set of $\{{\rho(t)}\}$ to extract an “experimental MLP” which may be compared directly with theory. These go as follows:

We begin by imposing the final boundary condition, i.e. we post–select on the desired $\rho_{T}$, at a later time $T$. This means that we must pick a distance measure $D(\rho_{1},\rho_{2})$ between quantum states (e.g. fidelity, Bures distance, or similar), and keep the sub–ensemble for which $D(\rho(t=T),\rho_{T})\leq W$, where $W$ is some small widow or allowed tolerance about the chosen final state. We’ll call the post–selected sub–ensemble (PSSE) $\{\rho(t)\}_{ps}$; this is, by construction, the set of trajectories which connect $\rho_{0}$ to $\rho_{T}$. 2. 2.

A simple intuition about the meaning of the MLP is that it should follow a densest cluster of trajectories in $\{\rho(t)\}_{ps}$; in order to approximate this concept of a “densest cluster” numerically, we must rank each SQT in $\{\rho(t)\}_{ps}$ according to its distance to all other SQTs in $\{\rho(t)\}_{ps}$. It is useful to construct a matrix of elements $\mathfrak{D}_{nm}$, where $n$ and $m$ are indices which run over the SQTs in $\{\rho(t)\}_{ps}$; we write

[TABLE]

where the sum over $k$ runs over all the timesteps between $t=0$ and $t=T$. The matrix $\mathfrak{D}$ will be symmetric as long as the distance measure $D$ is symmetric (we strongly discourage the use of any asymmetric distance measures for our present purposes). Then each SQT can be assigned a distance score relative to all other elements of $\{\rho(t)\}_{ps}$ according to $\bar{\mathfrak{D}}_{n}=\sum_{m}\mathfrak{D}_{nm}$, where a relatively smaller value of $\bar{\mathfrak{D}}_{n}$ indicates that trajectory $n$ of is closer to other trajectories in the post–selected set. These distances scores thus allow us to rank all of the trajectories in the PSSE. 3. 3.

The final step in the procedure is a simple average; we take the the closest–clustered 5%-10% of trajectories in $\{\rho(t)\}_{ps}$, (those with the smallest 5%-10% of $\bar{\mathfrak{D}}_{n}$), and average them. The idea is that this approximates the smooth curve following the densest cluster of SQTs in $\{\rho(t)\}_{ps}$.

We apply this procedure, and compare with the analytic solutions to (70), for a few selected boundary conditions, in Fig. 8.

We close with a few remarks about the procedure we have just described. We lack a formal proof that the numerical procedure just outlined necessarily always converges to the optimization we perform by the CDJ path integral method Chantasri et al. (2013). We nonetheless see in Fig. 8 that the present case continues to support the agreement between SQTs (from either simulation, or real data) and OPs, which have been successfully compared in many other scenarios as well Weber et al. (2014); Naghiloo et al. (2017); Lewalle et al. (2018). An important feature of the numerical procedure we have described is that it ranks trajectories by considering their entire evolution connecting $\mathbf{q}_{i}$ and $\mathbf{q}_{f}\tilde{\pm}W$, rather than looking for some piecewise optimization, or explicitly following a peak trajectory density at particular times. This reflects the character of the variational approach we have used in the theory; invoking $\delta\mathcal{S}=0$ leads to an optimization over an entire trajectory duration. While in many simple cases some piecewise numerical optimization could get us similar results, we would expect this simplified picture to fail in other cases, as it makes a substantial conceptual departure from the theory we want these numerics to match. In terms of practical considerations, we point out that the “experimental MLP” procedure above leads to an attrition in the number of SQTs used at each step, first due to post–selection, and then due to the ranking procedure. In order to obtain smooth results, it is necessary that at least several hundred SQTs make it into the final, most–closely clustered group which are averaged in step 3 above. The number of trajectories before post–selection, the size of the post–selection window, and the overall probability for trajectories to reach the desired final state, all play a role in these final numbers. The finer the post–selection window, the more initial trajectories are required at the beginning in order to obtain smooth results. Post–selections on very rare events may be prohibitively difficult to verify in practice, simply due to the overwhelming amount of data that would need to be collected / generated in order to have an appropriate number of SQTs in the PSSE.