Linear feedback stabilization of a dispersively monitored qubit

TL;DR

This paper demonstrates a feedback control method that stabilizes a dispersively monitored qubit near a target state by linearly feeding measurement records into the drive, accounting for real-world experimental imperfections.

Contribution

It introduces a novel linear feedback scheme for dispersively monitored qubits that enhances state stabilization, including analysis of practical nonidealities.

Findings

Feedback effectively stabilizes the qubit state near the target.

Feedback delay significantly impacts stabilization performance.

Optimal stabilization occurs when collapse timescale exceeds feedback delay.

Abstract

The state of a continuously monitored qubit evolves stochastically, exhibiting competition between coherent Hamiltonian dynamics and diffusive partial collapse dynamics that follow the measurement record. We couple these distinct types of dynamics together by linearly feeding the collected record for dispersive energy measurements directly back into a coherent Rabi drive amplitude. Such feedback turns the competition cooperative, and effectively stabilizes the qubit state near a target state. We derive the conditions for obtaining such dispersive state stabilization and verify the stabilization conditions numerically. We include common experimental nonidealities, such as energy decay, environmental dephasing, detector efficiency, and feedback delay, and show that the feedback delay has the most significant negative effect on the feedback protocol. Setting the measurement collapse timescale to be long compared to the feedback delay yields the best stabilization.

Figures

Equations

$${\hat{H}}_c=\hslash \left[{\Delta }_0+{\Delta }_1\tilde{r}(t-T_d)\right]\frac{{\hat{\sigma }}_{\varphi }}{2},$$

$${\hat{M}}_{\bar{r}}=\sqrt{P(\bar{r}|0)}|0⟩⟨0|+\sqrt{P(\bar{r}|1)}|1⟩⟨1|.$$

$${\hat{M}}_{\bar{r}}$$

$$\hat{U}$$

$$\hat{\rho }(t+dt)=({\mathcal{E}}_{\gamma }\circ {\mathcal{E}}_{T_2}\circ {\mathcal{E}}_{T_1})\left[\frac{\hat{U}{\hat{M}}_{\bar{r}}\hat{\rho }(t){\hat{M}}_{\bar{r}}^{†}{\hat{U}}^{†}}{\mathrm{Tr}\left({{\hat{M}}_{\bar{r}}^{†}\hat{M}}_{\bar{r}}\hat{\rho }(t)\right)}\right],$$

$$P(\bar{r})=\text{Tr}\left({\hat{M}}_{\bar{r}}^{†}{\hat{M}}_{\bar{r}}\hat{\rho }\right)$$

$$\approx \frac{\exp (-dt(\bar{r}-z{)}^2/2{\tau }_m)}{\sqrt{2\pi {\tau }_m/dt}}.$$

$$p_n$$

$$x_n^{\prime }$$

$$y_n^{\prime }$$

$$z_n^{\prime }$$

$$\Delta$$

$$y_n^{\prime \prime }$$

$$z_n^{\prime \prime }$$

$$x_{n+1}$$

$$y_{n+1}$$

$$z_{n+1}$$

$$\bar{r}(t)=z(t)+\sqrt{{\tau }_m}\xi (t)$$

$$d\hat{\rho }$$

$$+\frac{{\Gamma }^{\prime }}{2}\mathcal{D}[{\hat{\sigma }}_z]\hat{\rho }dt+\frac{1}{T_1}\mathcal{D}[{\hat{\sigma }}_{-}]\hat{\rho }dt+\frac{{\tau }_m{\Delta }_1^2}{4}\mathcal{D}[{\hat{\sigma }}_x]\hat{\rho }dt$$

$$+\frac{1}{2}\mathcal{H}[{\hat{\sigma }}_z]\hat{\rho }\frac{dW}{\sqrt{{\tau }_m}}+i\frac{{\tau }_m{\Delta }_1}{2}[{\hat{\sigma }}_x,\hat{\rho }]\frac{dW}{\sqrt{{\tau }_m}}.$$

$$\dot{x}$$

$$\dot{y}$$

$$-yz\frac{\xi }{\sqrt{{\tau }_m}}+{\tau }_m{\Delta }_{1}z\frac{\xi }{\sqrt{{\tau }_m}},$$

$$\dot{z}$$

$$+(1-z^2)\frac{\xi }{\sqrt{{\tau }_m}}-{\tau }_m{\Delta }_{1}y\frac{\xi }{\sqrt{{\tau }_m}},$$

$${\Delta }_0$$

$$F=⟨{\theta }_s|\hat{\rho }|{\theta }_s⟩=\frac{1+R}{2}$$

$$y_s$$

$$z_s$$

$$\tan {\theta }_s$$

$$\stackrel{T_1\to \infty }{\to }-\frac{{\tau }_m{\Delta }_1^2}{2{\Delta }_0},$$

$$R_s$$

$$\stackrel{T_1\to \infty }{\to }\frac{|{\Delta }_1|}{|{\Delta }_0||\sec {\theta }_s|+\Gamma |\sin {\theta }_s|}$$

$${\Delta }_0$$

$${\Delta }_1$$

$$R_{\max }(\theta )$$

$$\stackrel{T_1\to \infty }{\to }\frac{1}{\sqrt{2{\tau }_m\Gamma }}=\frac{1}{\sqrt{\frac{1}{\eta }+\frac{2{\tau }_m}{T_2}}}.$$

$${\Delta }_0$$

$$\stackrel{T_1\to \infty }{\to }-\frac{\sin 2\theta }{4{\tau }_m}\left[\frac{1}{\eta }+\frac{2{\tau }_m}{T_2}\right]$$

$${\Delta }_1$$

$${\widetilde{\stackrel{ˉ}{r}}}_k$$

Full text

I Introduction

Modern efforts to build a quantum computer have recently enabled the time-continuous measurement of quantum state trajectories Gardiner (2004); Carmichael (1993); Wiseman and Milburn (2010); Jacobs and Steck (2006) using superconducting circuit quantum-electrodynamics (cQED) Gambetta et al. (2008); Korotkov (2014, 2016). For quantum non-demolition (QND) measurements, such as dispersive transmon energy level measurements Koch et al. (2007), a time-continuous noisy signal gradually reveals information over time about the stationary eigenstates of the measurement. Learning this information causes a gradual collapse of the qubit state over time to one of the eigenstates of the measurement, which can be observed Murch et al. (2013); Roch et al. (2014); Campagne-Ibarcq et al. (2016a) by processing the noisy measurement record. Notably, such gradual collapse is not monotonic, but rather stochastic, with some temporal increments effectively erasing the information gathered by previous increments and “uncollapsing” the state Korotkov and Jordan (2006); Katz et al. (2008); Jordan and Korotkov (2010); Schindler et al. (2013). Nevertheless, the random walk of the qubit state does eventually fully collapse to one of the measurement eigenstates, where it remains. Continuing to monitor an eigenstate after collapse helps protect it from environmental decay due to the quantum Zeno effect Misra and Sudarshan (1977); Facchi and Pascazio (2008); Slichter et al. (2016), which suppresses other coherent evolution via continual recollapse.

The stochastic collapse evolution from continuous monitoring generally competes with other coherent Hamiltonian evolution, which produces complicated state dynamics that have rich statistical correlations Weber et al. (2014); Chantasri et al. (2013); Jordan et al. (2016); Chantasri and Jordan (2015); Naghiloo et al. (2017); Jordan and Büttiker (2005); Chantasri et al. (2016); Lewalle et al. (2017). Similarly, simultaneously monitoring multiple non-commuting observables produces nontrivial competitive dynamics Ruskov et al. (2010); Hacohen-Gourgy et al. (2016); García-Pintos and Dressel (2016); Atalaya et al. (2017). In either case, collapse to a stable eigenstate is prevented. Remarkably, it has been shown that such dynamical competition can be made cooperative by feeding a suitably processed noisy readout back into the controlled Hamiltonian dynamics Wiseman and Milburn (1993); Wiseman (1994); Hofmann et al. (1998); Doherty and Jacobs (1999); Wang and Wiseman (2001); Wiseman et al. (2002); Ruskov and Korotkov (2002); Ruskov et al. (2004); Zhang et al. (2005); Korotkov (2005); Wang et al. (2005); Schirmer and Wang (2010); Chia and Wiseman (2011). Similar to classical control theory, the random perturbations induced by the stochastic dynamics can be compensated by the adaptive control Hamiltonian to produce customized dynamics that approximate a desired outcome. Unlike classical control theory, such measurement-based quantum feedback control can fundamentally leverage the peculiarities of quantum state collapse in ways that have no classical analog Doherty et al. (2000). For example, quantum states may be more rapidly purified through a clever application of feedback Jacobs (2003); Combes and Jacobs (2006); Wiseman and Ralph (2006); Jordan and Korotkov (2006); Combes et al. (2008); Ruskov et al. (2012). A variety of quantum feedback control schemes have been implemented within the past six years, including photon number state preparation Sayrin et al. (2011), continuous superconducting Rabi oscillation stabilization Vijay et al. (2012), discrete feedback control of superconducting qubits Risté et al. (2012); Campagne-Ibarcq et al. (2013), entangled state generation Risté et al. (2013); Shankar et al. (2013); Lin et al. (2013); Leghtas et al. (2015); Schwartz et al. (2016), discrete reversal of continuous random phase drift de Lange et al. (2014), and continuous qubit state stabilization via both dissipative bath engineering Murch et al. (2012) and feedback from fluorescence measurements Campagne-Ibarcq et al. (2016b).

In this paper, we theoretically revisit the latter task of state stabilization, but for a transmon qubit undergoing continuous linear feedback of a dispersive energy measurement record in a realistically modeled environment that includes dephasing, energy-decay, measurement inefficiency, signal filtering, and feedback delay. The goal of such a stabilization protocol is to effectively alter the eigenstates of the measurement collapse to target an arbitrary desired state, while preserving the main features of the quantum Zeno effect that protect that state against deterioration from environmental factors. Unlike previously considered dissipative bath engineering Murch et al. (2012) or fluorescence-based feedback Wang and Wiseman (2001); Campagne-Ibarcq et al. (2016b) methods, such dispersive feedback is based on a fundamentally QND measurement. As such, the symmetry of the non-feedback collapse dynamics alters the permissible state stabilization regions dramatically from those produced by the asymmetric backaction of environmental dissipation or energy-decay from fluorescence. The QND nature of the dispersive measurement also makes it natural to compare the differences between ensemble-average state stabilization and individual trajectory stabilization. We thus pose the stabilization problem in two contrasting ways: First, we suppose an experimenter only knows the applied feedback parameters and does not monitor individual measurement records, and so can only prepare a targeted ensemble-averaged state on demand using feedback. Second, we suppose the experimenter also processes each measurement record to have more detailed information about individual quantum state trajectories and so can choose whether a particular state is suitably prepared on demand. We include environmental effects in our analysis with parameters chosen to be consistent with current experimental technology for dispersive transmon trajectory measurements Murch et al. (2013); Weber et al. (2014); Roch et al. (2014); Hacohen-Gourgy et al. (2016).

In the Markovian case with no signal filtering or feedback delay, we solve the two variations of the problem analytically and verify the results numerically. We find that in the ideal case with no environmental effects, both ways of posing the problem coincide: Any angle in the qubit Bloch sphere may be perfectly targeted by a suitable choice of linear feedback parameters. When environmental effects are added (neglecting feedback delay and signal filtering), we find that the two variations still coincide provided the feedback parameters are optimally chosen. That is, targeting any angle in the Bloch sphere with maximum ensemble-averaged purity is equivalent to stabilizing individual trajectories such that their measurement backaction is minimized. We find that all angles not too near the natural measurement pole corresponding to the excited qubit state may still be targeted with high purity, and solve for the feedback parameters needed to obtain each such state. Notably, while the ensemble-averaged purity remains nearly constant, the distributions of individual trajectories spread asymmetrically near the measurement poles, developing peaks of maximum probability with higher purity than the ensemble average. For stabilized states at the equator of the Bloch sphere, the maximum probability peak coincides with the ensemble-averaged state. For stabilized states near the excited state measurement pole, the maximum probability peak bifurcates into two distinct regions, with one dominant, such that each region is substantially more pure than the ensemble-average state. As a result, the two ways of posing the stabilization problem have distinct solutions for states near the natural measurement poles. For stabilized states near the ground state measurement pole, both the ensemble average and peak purify due to energy-relaxation environmental effects.

Once signal filtering and feedback delay are included, the dynamics become intrinsically non-Markovian due to the delay buffer. We investigate these effects numerically. Compared to the Markovian case, we find that the stabilization is minimally degraded by single-pole (RC) low-pass filtering of the signal, so the protocol appears robust to realistic frequency-filtering effects in the feedback circuitry. We also find that feedback delay has the largest negative effect on the protocol, causing a dramatic reduction in stabilization purity as the delay becomes comparable to the measurement collapse timescale. This sensitivity to delay agrees with the feedback analysis performed for stabilized Rabi oscillations of a double-quantum-dot setup in Ref. Ruskov et al. (2004). In our case, the trajectory distributions rapidly broaden with delay since the delayed signal being fed back into the controller has diminished relevance for the later state evolution that is being controlled. Nevertheless, the effect of delay can be mitigated by slowing the measurement rate to keep the collapse timescale an order of magnitude longer than the delay in the feedback circuitry. Curiously, despite the fact that in the Markovian case the Bloch sphere equator can be stabilized, we find that any amount of filtering or delay destabilizes a sharp angular region around the equator, which then becomes consistent with the instability observed at the equator in fluorescence-based protocols Wang and Wiseman (2001); Campagne-Ibarcq et al. (2016b). In summary, most angles in the qubit Bloch sphere may be stabilized with high purity even when realistic models of the environment and feedback circuitry are taken into account.

The paper is organized as follows. In Sec. II, we detail the modeling of the dispersive qubit monitoring with feedback for both idealized and realistic situations. In Sec. III, we pose the problem of ensemble-averaged qubit state stabilization and discuss both analytical results and numerical simulations. In Sec. IV, we pose the related problem of state stabilization for individual qubit trajectories and discuss both analytical results and numerical simulations. We conclude in Sec. V.

II Feedback Model

For concreteness, we consider a superconducting transmon qubit that is dispersively coupled to a microwave resonator as the measuring apparatus Koch et al. (2007); Gambetta et al. (2008); Korotkov (2014, 2016). The qubit Hamiltonian describing the lowest two levels of the transmon is $\hat{H}_{q}=\hbar\omega_{q}\,\hat{\sigma}_{z}/2$, where the Pauli operator $\hat{\sigma}_{z}=\ket{1}\!\bra{1}-\ket{0}\!\bra{0}$ indicates the difference between the excited ($\ket{1}$) and ground ($\ket{0}$) state energy levels, and $\omega_{q}$ is the natural qubit frequency. We will work in the rotating frame of the qubit in what follows for simplicity to neglect this natural evolution. The Hamiltonian of the coupled microwave resonator mode is harmonic $\hat{H}_{r}=\hbar\omega_{r}\,\hat{a}^{\dagger}\hat{a}$, where $\hat{a}$ is its lowering operator and $\omega_{r}$ is the natural resonator frequency that is detuned from the qubit. The dispersive coupling Hamiltonian between the qubit and resonator then approximates $\hat{H}_{qr}=\hbar\chi\,\hat{\sigma}_{z}\,\hat{a}^{\dagger}\hat{a}$, where $\pm\chi$ is the qubit-state-dependent frequency shift of the resonator due to the dispersive coupling. The resonator has energy-decay rate $\kappa$ and is pumped on resonance with a coherent field to reach a qubit-state-dependent coherent steady state with qubit-state-independent average photon number $\bar{n}$. The steady-state microwave field then leaks from the resonator and travels down a transmission line, where it is amplified along the maximally informative quadrature Korotkov (2014, 2016) with a phase-sensitive amplifier and finally measured to produce a time-continuous homodyne signal $I(t)$. Due to the entanglement with the leaked microwave field, definite qubit states $\ket{0}$ and $\ket{1}$ produce noisy homodyne signals that temporally average to distinct mean values $I_{0}$ and $I_{1}$, respectively. As such, monitoring the homodyne signal provides time-continuous information about the qubit energy basis. On average, these homodyne measurements decohere the qubit state at the ensemble-averaged measurement-dephasing rate $\Gamma_{m}=8\chi^{2}\bar{n}/\kappa$ Gambetta et al. (2008); Korotkov (2014, 2016), which may be determined experimentally by comparing the Ramsey decay rates of the qubit with and without measurement Murch et al. (2013); Weber et al. (2014); Roch et al. (2014); Chantasri et al. (2016).

To stabilize the qubit state, we process the homodyne signal $I(t)$ and feed it back into the coherent qubit evolution as a coherent Rabi drive Wiseman and Milburn (1993); Wiseman (1994); Hofmann et al. (1998); Doherty and Jacobs (1999); Wang and Wiseman (2001); Wiseman et al. (2002); Wang et al. (2005); Campagne-Ibarcq et al. (2016b). For convenience, we first rescale the signal so that it directly corresponds to $\hat{\sigma}_{z}$. The rescaled readout has the form $r(t)\equiv 2(I(t)-\bar{I})/\Delta I$, where $\bar{I}=(I_{0}+I_{1})/2$ is the mean homodyne current, and $\Delta I=I_{1}-I_{0}$ is the homodyne signal contrast. The rescaled readout $r(t)$ then has a temporal average of $\pm 1$ for definite qubit energy states, in agreement with the eigenvalues of $\hat{\sigma}_{z}$. Before being fed back into the controller, this rescaled readout is necessarily frequency-filtered by the feedback circuitry Ruskov et al. (2004), $r\mapsto\tilde{r}$, which also delays the filtered readout by a duration $T_{d}$. This filtered and delayed readout $\tilde{r}(t-T_{d})$ is then used to modulate a coherent Rabi drive described by the control Hamiltonian,

[TABLE]

where $\hat{\sigma}_{\phi}=\cos\phi\,\hat{\sigma}_{x}+\sin\phi\,\hat{\sigma}_{y}$, with $\hat{\sigma}_{x}=\ket{1}\!\bra{0}+\ket{0}\!\bra{1}$ and $\hat{\sigma}_{y}=-i\ket{1}\!\bra{0}+i\ket{0}\!\bra{1}$. The choice of angle $\phi$ fixes the plane of the qubit Bloch sphere to which the qubit will become stabilized. For simplicity of discussion, we now fix $\phi=\pi$ to choose clockwise oscillations in the $yz$-plane. This feedback Hamiltonian induces oscillations at a constant Rabi frequency $\Delta_{0}$ that is modulated by the linear feedback term $\Delta_{1}\,\tilde{r}(t-T_{d})$. As we will show, and in agreement with fluorescence-based feedback stabilization protocols Wang and Wiseman (2001); Wiseman et al. (2002); Wang et al. (2005); Campagne-Ibarcq et al. (2016b), the two parameters $\Delta_{0}$ and $\Delta_{1}$ completely control the effect of the feedback by jointly determining the attraction region in state space for stabilization. In contrast to the fluorescence-based protocols that have a bias toward the ground state, we will show that dispersive measurements can stabilize nearly any angle within the chosen plane with high purity (i.e., nearly any pure state may be targeted with high fidelity).

We now describe the qubit evolution produced by collecting the dispersive readout $r(t)$ and applying the feedback Hamiltonian in Eq. (1). We first consider a time-discrete quantum Bayesian model Korotkov (2014, 2016), from which we can numerically simulate the feedback protocol including inefficiencies and delay. We then consider a time-continuous stochastic master equation model Gambetta et al. (2008) that interpolates the Markovian feedback dynamics in the limit of zero delay $T_{d}\to 0$ (and no signal filtering) with a fictitious time-continuous stochastic process Wiseman and Milburn (1993); Wiseman (1994); Hofmann et al. (1998); Wang and Wiseman (2001); Wiseman et al. (2002); Wang et al. (2005), and which is convenient for analytic derivations of the required feedback parameters.

II.1 Time-discrete model

The strength of the qubit measurement is determined by the timescale $\tau_{m}$ needed to obtain a unit signal-to-noise ratio for distinguishing the qubit states. A qubit initially in a superposition of energy states will typically collapse to a definite energy state within a few $\tau_{m}$. This collapse timescale $\tau_{m}$ is related to the ensemble dephasing rate $\Gamma_{m}$ of the qubit due to measurement according to $\Gamma_{m}=(2\tau_{m}\eta)^{-1}$, where $\eta\in[0,1]$ is the quantum efficiency of the microwave detection circuitry Korotkov (2014, 2016).

When the energy decay rate of the resonator $\kappa$ is sufficiently fast compared to both the measurement rate $\Gamma_{m}$ and the dynamics of the qubit (known as the “bad cavity regime”), the coherently pumped resonator effectively remains at steady-state. We may then approximately ignore the resonator and phenomenologically describe the qubit evolution alone Korotkov (2014), provided that we coarse-grain the collected readout $r(t)$ into discrete segments $\{\bar{r}_{k}\}$ that are temporal averages $\bar{r}_{k}\,dt\equiv\int_{t_{k}}^{t_{k}+dt}r(t^{\prime})dt^{\prime}$ over time increments $dt>\kappa^{-1}$ longer than the relaxation timescale of the resonator field Korotkov (2016). This coarse-grained measurement evolution is then Markovian in the discretized time intervals, and can be efficiently described using the quantum Bayesian method Korotkov (2014, 2016), which also permits efficient numerical simulation.

The total evolution of the qubit can be approximately decomposed into three parts: pure measurement backaction, coherent unitary evolution, and environmental dissipation. These three contributions may be treated separately provided that the duration $dt$ is smaller than the characteristic timescales of each evolution. We now describe each contributing evolution, then combine them.

First, the form of the measurement backaction follows from Bayes’ theorem and the Born rule Korotkov (2014). Given empirically measured conditional probabilities $P(\bar{r}|0)$ and $P(\bar{r}|1)$ for obtaining a result $\bar{r}$ averaged over $dt$ given a definite qubit energy state (0 or 1), the qubit energy probabilities must update according to Bayes’ theorem $P(0)\mapsto P(0|\bar{r})=P(\bar{r}|0)P(0)/P(\bar{r})$ and $P(1)\mapsto P(1|\bar{r})=P(\bar{r}|1)P(1)/P(\bar{r})$, where the total probability for collecting $\bar{r}$ is $P(\bar{r})=P(\bar{r}|0)P(0)+P(\bar{r}|1)P(1)$. This evolution is equivalent to the density operator update rule $\hat{\rho}\mapsto\hat{M}_{\bar{r}}\hat{\rho}\hat{M}^{\dagger}_{\bar{r}}/P(\bar{r})$ with probability $P(\bar{r})=\text{Tr}(\hat{M}^{\dagger}_{\bar{r}}\hat{M}_{\bar{r}}\hat{\rho})$, determined entirely by the Kraus operator

[TABLE]

For a maximally informative choice of amplified quadrature, there is no additional phase-backaction on the qubit Korotkov (2014, 2016), so this minimal Kraus operator is sufficient.

Because of the central limit theorem, the conditional probabilities $P(\bar{r}|0)$ and $P(\bar{r}|1)$ approximate Gaussians centered at $\pm 1$. For simplicity, we also assume approximately equal variances, which may be justified empirically. To preserve proper scaling of temporal averages as the discrete time step $dt$ is varied, the variances of each time step must be inversely proportional to $dt$. By definition of the unit signal-to-noise timescale $\tau_{m}$, the variances must then have the form $\tau_{m}/dt$, yielding

[TABLE]

where $C(\bar{r})=(dt/2\pi\tau_{m})^{1/4}e^{-(\bar{r}^{2}+1)dt/4\tau_{m}}$ is a state-independent normalization factor that cancels in the state update rule. Note that this Kraus operator may be conveniently written $\hat{M}_{\bar{r}}=C(\bar{r})\exp(\bar{r}dt\hat{\sigma}_{z}/2\tau_{m})=C(\bar{r})[\cosh(\bar{r}dt/2\tau_{m})\hat{1}+\sinh(\bar{r}dt/2\tau_{m})\hat{\sigma}_{z}]$, with $\hat{1}=\ket{0}\!\bra{0}+\ket{1}\!\bra{1}$. The rate $(2\tau_{m})^{-1}$ that appears here is the ensemble-average dephasing rate obtained after averaging over all collected signal $\bar{r}$.

Second, the form of the unitary dynamics follows from the control Hamiltonian $\hat{H}_{c}$ in Eq. (1) with $\phi=\pi$, which describes precession in the $yz$-plane of the qubit Bloch sphere. Treating the coherent evolution independently yields the unitary operator $\hat{U}=\exp(-idt\hat{H}_{c}/\hbar)$ that simplifies to the convenient form

[TABLE]

with a frequency $\Delta=\Delta_{0}+\Delta_{1}\,\tilde{\bar{r}}(t-T_{d})$ that depends upon the filtered average readout $\tilde{\bar{r}}(t-T_{d})$ collected from a time delay $T_{d}=n_{d}\,dt$ of $n_{d}$ time steps in the past and fed back into the controller during the current time interval $dt$. The nonzero feedback delay makes the resulting dynamics intrinsically non-Markovian. Additional frequency-filtering from the feedback circuitry Ruskov et al. (2004) is discretely modeled by a recursive transformation of the history of raw average readouts $\{\bar{r}_{k}\}$ stored in a memory buffer (see Sec. III.3 for more detail).

Third, we consider the most common forms of environmental dissipation: energy relaxation, energy dephasing, and detector inefficiency. Energy relaxation follows a phenomenological decay rate $1/T_{1}$ for the excited state to relax into the ground state. Defining the density operator matrix elements as $P_{0}=\bra{0}\hat{\rho}\ket{0}$, $P_{1}=\bra{1}\hat{\rho}\ket{1}$, and $\rho_{01}=\bra{0}\hat{\rho}\ket{1}=\rho_{10}^{*}\propto\sqrt{P_{0}P_{1}}$, energy decay yields the transformation $P_{1}\mapsto P_{1}\exp(-dt/T_{1})$, which forces $P_{0}\mapsto 1-P_{1}\exp(-dt/T_{1})$ and $\rho_{01}\mapsto\rho_{01}\exp(-dt/2T_{1})$. We summarize these rules with the energy-decay map $\hat{\rho}\mapsto\mathcal{E}_{T_{1}}(\hat{\rho})$. Similarly, environmental energy dephasing follows a phenomenological decay rate $1/T_{2}$, and corresponds to the transformation $\rho_{01}\mapsto\rho_{01}\exp(-dt/T_{2})$, summarized by the dephasing map $\hat{\rho}\mapsto\mathcal{E}_{T_{2}}(\hat{\rho})$. The detector inefficiency $(1-\eta)$ describes the fraction of measured qubit information that has been lost to the environment, thus producing residual measurement dephasing on average with rate $\gamma=\Gamma_{m}-1/2\tau_{m}=(1-\eta)/(2\tau_{m}\eta)$, which is summarized by the map $\hat{\rho}\mapsto\mathcal{E}_{\gamma}(\hat{\rho})$ that dephases in a similar way to $\mathcal{E}_{T_{2}}$; i.e., $\rho_{01}\mapsto\rho_{01}\exp(-\gamma dt)$.

The total state update is then approximately described by the composite map

[TABLE]

where all operations are described at time $t$. This separation of the evolution into distinct pieces is formally valid to linear order in $dt$, but (unlike explicitly linear-order updates) ensures completely positive evolution of the state.

To minimize accumulated error from the composition approximation in Eqs. (23) and (5), the time step should satisfy $dt/\tau_{m}\ll 1$, $dt\,\Delta\ll 1$, and $dt\,\Gamma_{m}\ll 1$. In particular, since the noisy feedback signal $\tilde{\bar{r}}(t-T_{d})$ will likely dominate the rate $\Delta=\Delta_{0}+\Delta_{1}\tilde{\bar{r}}(t-T_{d})$, it should satisfy $dt\,\tilde{\bar{r}}\,\Delta_{1}\ll 1$. Since $\bar{r}\lesssim 5\,\sqrt{\tau_{m}/dt}$ by Gaussian statistics, this implies $dt\,\tilde{\bar{r}}\,\Delta_{1}<5\Delta_{1}\sqrt{dt\,\tau_{m}}\ll 1$, so $\Delta_{1}\ll 1/(5\sqrt{dt\,\tau_{m}})$ is a practical upper bound for $\Delta_{1}$ for a chosen time step $dt$ in Eq. (5). We also note that the chosen ordering of operations, in particular measurement backaction followed by feedback control, is made to anticipate the Markovian feedback limit with vanishing delay ($T_{d}\to 0$) in the next section.

To efficiently simulate Eq. (5) numerically, we use a Bloch coordinate decomposition of the state, $\hat{\rho}\to(x,y,z)$, where $x=\text{Tr}\left(\hat{\sigma}_{x}\,\hat{\rho}\right)$, $y=\text{Tr}\left(\hat{\sigma}_{y}\,\hat{\rho}\right)$, and $z=\text{Tr}\left(\hat{\sigma}_{z}\,\hat{\rho}\right)$. With this representation, the rescaled average readout $\bar{r}$ at each time $t$ may be sampled from the readout distribution

[TABLE]

where the single-Gaussian approximation makes numerical sampling efficient. After generating each $\bar{r}$, the update in Eq. (5) is applied, using the previously sampled, time-delayed, and filtered $\tilde{\bar{r}}$ as appropriate in the feedback unitary.

Expressed in Bloch-coordinates, Eq. (5) has a simple prescription: Given the state $(x_{n},y_{n},z_{n})$ and a random $\bar{r}$ sampled at time step $t_{n}$, as well as a filtered $\tilde{\bar{r}}_{n-n_{d}}$ from time step $t_{n}-T_{d}$ with $T_{d}=n_{d}\,dt$, the new state $(x_{n+1},y_{n+1},z_{n+1})$ at time $t_{n+1}=t_{n}+dt$ is obtained by the following sequence of transformations

[TABLE]

Note that the normalization factor $p_{n}$ is proportional to the probability for sampling $\bar{r}$, but irrelevant constants have been canceled. At the start of each simulated evolution, there is no pre-sampled readout history, so we treat the applied feedback $\tilde{\bar{r}}_{n-n_{d}}$ as zero in the update until the feedback delay buffer has been filled. For simplicity, we set $x=0$ initially, which constrains the dynamics to the $yz$-plane.

II.2 Time-continuous model

For analytic convenience, the time-continuum limit $dt\to 0$ can be used to derive stochastic master equations (SMEs) from the update equation Eq. (5) in the Markovian feedback limit with $T_{d}\to 0$. Note that this limit formally produces a physically incorrect model for two reasons: First, $dt>\kappa^{-1}$ must be satisfied to phenomenologically ignore the relaxation of the resonator to its steady state. The continuum limit thus implies a resonator with infinitely fast energy decay, $\kappa\to\infty$, which would produce no qubit measurement, since $\Gamma_{m}\propto\kappa^{-1}\to 0$. Second, $T_{d}>0$ in any realistic laboratory feedback control loop due to finite signal velocity and controller processing time, so the limit $T_{d}\to 0$ is unrealizable. Nevertheless, this Markovian feedback continuum limit still formally produces a useful interpolation of the discretized dynamics in the preceding section as an idealized time-continuous stochastic process that may be more easily analyzed Gambetta et al. (2008). That is, averaging this fictitious process over time bins $dt>\kappa^{-1}$ should still recover the physically correct time-discretized evolution. In particular, it is simple to obtain a master equation (ME) describing the ensemble-averaged behavior from such an SME. It is also a good approximation when the delay is short compared to the intrinsic collapse time of the measurement, $T_{d}\ll\tau_{m}$. Indeed, such a separation of time scales is always implied by such time-continuous stochastic equations of motion Gardiner (2004); Carmichael (1993).

Stochastic processes are not differentiable. Thus, a careful treatment is needed to derive the differential equations for the evolution Gardiner (2004). There are two popular pictures for stochastic differential equations, Itô and Stratonovich, the former having the advantage that the stochastic noise at any time $t$ is independent from the process in which it occurs (see Ref. Gardiner (2004) and Appendix B of Ref. Wiseman and Milburn (2010)). This advantage makes finding the ME for ensemble averaged evolution straightforward in the Itô picture, so we opt for that approach here.

Following the Itô formalism, we observe that Eq. (II.1) may be written as $P(\bar{r})=P(dW)dt/\sqrt{\tau_{m}}$, with the zero-mean Gaussian distribution $P(dW)=\exp(-dW^{2}/2dt)/\sqrt{2\pi dt}$ having variance $dt$, provided that we make the identification $\bar{r}\,dt=z\,dt+\sqrt{\tau_{m}}\,dW$. The random variable $dW$ then has the properties of a Wiener increment, so that the identity $dW^{2}=dt$ becomes exact in the continuum limit $dt\to 0$ Gardiner (2004). As such, we can interpolate the coarse-grained time-dependent readout as a moving-mean stochastic process

[TABLE]

centered at $z(t)$ with additive white noise $\xi(t)=dW/dt$ with zero-mean $\langle\xi(t)\rangle=0$ and Dirac-$\delta$ temporal-correlation $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$.

We can now expand the right hand side of the update equation Eq. (5) to first order of $dt$ while adhering to the Itô rule $dW^{2}=dt$, and assuming $T_{d}\to 0$ (implying the filtered readout $\tilde{\bar{r}}$ in $\hat{U}$ becomes identified with the just-sampled and unfiltered $\bar{r}$ in $\hat{M}_{\bar{r}}$). This expansion produces the Markovian feedback Itô stochastic master equation

[TABLE]

Here $\Gamma^{\prime}=1/T_{2}+1/2\tau_{m}\eta$ and $\hat{\sigma}_{-}=\ket{0}\!\bra{1}$. We also define the dissipation superoperator $\mathcal{D}[\hat{A}]\hat{\rho}=\hat{A}\hat{\rho}\hat{A}^{\dagger}-\{\hat{A}^{\dagger}\hat{A},\,\hat{\rho}\}/2$, as well as the innovation superoperator $\mathcal{H}[\hat{A}]\hat{\rho}=\hat{A}\hat{\rho}+\hat{\rho}\hat{A}^{\dagger}-\text{Tr}(\hat{A}\hat{\rho}+\hat{\rho}\hat{A}^{\dagger})\hat{\rho}$ Wiseman and Milburn (2010). Of particular importance is the commutator of the anti-commutator that appears in the first line, which describes the effect of the feedback control applied immediately after the collapse induced by the measurement Wiseman and Milburn (1993); Wiseman (1994); Wiseman and Milburn (2010). The ensemble-averaged ME is obtained simply by setting the terms proportional to $dW$ to zero in Eq. (II.2).

This result simplifies when written in terms of the Bloch coordinates, yielding

[TABLE]

where $\Gamma=1/2T_{1}+1/T_{2}+1/2\tau_{m}\eta$ is the total ensemble $z$-dephasing rate in the absence of feedback. Note that the feedback contributes additional $x$-dephasing at a rate $\tau_{m}\Delta_{1}^{2}/2$. The ensemble-averaged ME can be recovered from this SME by setting the terms proportional to $\xi$ to zero. The consistency between the deterministic ensemble-average ME and the discrete formulation of Eq. (5) is then checked by averaging several thousand trajectories numerically generated using Eq. (5) and comparing to the analytic solution of the deterministic part of Eq. (10). We will primarily use the Bloch coordinate representation in what follows for analytical calculations.

III Average State Stabilization

We now pose the following state stabilization problem: Consider a closed, continuous measurement feedback loop that is to be used as a storage mechanism for preparing a target qubit state on demand. An experimenter will press a button at a unknown future time, terminating the feedback loop, in order to receive the prepared qubit state for immediate experimentation. Can the feedback parameters $\Delta_{0}$ and $\Delta_{1}$ be fixed such that a particular qubit state is reliably prepared on average?

We show that in the case of idealized Markovian feedback an arbitrary pure state may be prepared through feedback alone. We then generalize the protocol by adding realistic experimental nonidealities, including energy decay, environmental dephasing, and measurement inefficiency, showing that the best prepared states have degraded purity compared to the ideal case, but may still be meaningfully stabilized by the feedback loop. Finally, we discuss the significant effect of realistic signal filtering and feedback delay on the stabilization protocol, and show that good stabilization is still possible when the delay is sufficiently short compared to the measurement collapse timescale.

III.1 Stationary states for ideal Markovian feedback

For the case of a pure state ($T_{1},T_{2}\to\infty,\,\eta=1$) with no feedback delay ($T_{d}\to 0$), we use the Bloch Eqs. (10) to answer the question posed above. Specifically, consider the case where $x=0$ initially, so the qubit state lies in the $yz$-plane at $(0,y_{s},z_{s})$, with $y_{s}=\sin\theta_{s}$ and $z_{s}=\cos\theta_{s}$ determined entirely by the polar angle in the Bloch sphere. We then demand the condition that this initial state be a fixed point of the ensemble-averaged dynamics, meaning that $\dot{y}=\dot{z}=0$ in Eq. (10) with $\xi\to 0$. Solving this constraint yields the feedback parameters

[TABLE]

That is, for any choice of state angle $\theta_{s}$, there exist unique feedback parameters $\Delta_{0}$ and $\Delta_{1}$ that stabilize that state on average. These parameters are shown as functions of the stabilized state angle $\theta_{s}$ in Figure 1. The dependence of the linear feedback parameter $\Delta_{1}$ on $\sin\theta_{s}$ is sensible, since for $\theta_{s}=\{0,\pi\}$ the natural measurement poles are already fixed points for collapse. Similarly, it is sensible for the constant rate $\Delta_{0}$ to have dependence on $\sin 2\theta_{s}$, since the constant rotation biases the fixed point asymmetrically and should vanish for the backaction-symmetric stabilization points on the equator with $\theta_{s}=\pm\pi/2$.

Numerical simulations are in excellent agreement with this result for the ideal stationary state, with the feedback consistently isolating a unique stabilized state. The only exception to the uniqueness of the stabilization is that the active control vanishes at the natural measurement poles, since $\Delta_{0},\,\Delta_{1}\rightarrow 0$ as $\theta_{s}\rightarrow 0,\,\pi$. The absence of feedback at the poles permits stochastic collapse to either $\hat{\sigma}_{z}$-eigenstate. This limitation can be readily ameliorated in practice, however, by stabilizing a state near the pole of choice before removing feedback and permitting a final collapse step with high probability to the desired pole, which yields a protocol analogous to standard state heralding.

III.2 Stationary states for Markovian feedback

For realistic laboratory situations, we must add the effects of the energy decay time $T_{1}$, dephasing time $T_{2}$, and imperfect measurement efficiency $\eta$. These non-idealities decrease the state purity, so states in the $yz$-plane will be characterized by both a polar angle $\theta$ and a radius $R$, related to the purity by $\text{Tr}(\hat{\rho}^{2})=(1+R^{2})/2$, such that $y=R\sin\theta$ and $z=R\cos\theta$. Since a pure state is characterized entirely by the polar angle $\theta$, we may also consider the radius $R$ to be a measure of the average state fidelity

[TABLE]

between a target pure state and a mixed state at the same stabilized angle $\theta_{s}$. We will see that the stabilized angle $\theta_{s}$ may be reliably targeted, which makes $R$ a sensible measure of stabilization quality (as an aside, it is equivalent to the “synchronization degree” $D$ used in Ref. Zhang et al. (2005)).

As before, we demand that $\dot{y}=\dot{z}=0$ in Eq. (10) for the ensemble average when $\xi\to 0$. This stability condition produces the following set of stationary states that correspond to each particular choice of feedback parameters $\Delta_{0}$ and $\Delta_{1}$

[TABLE]

However, not every state in the $yz$-plane may be stabilized due to constraints on the average purity. This fact becomes more clear when we express the stationary states in polar form with $\tan\theta_{s}=y_{s}/z_{s}$ and $R_{s}^{2}=y_{s}^{2}+z_{s}^{2}$,

[TABLE]

The negligible energy-decay limit $T_{1}\to\infty$ is simpler to analyze, yet still corresponds to a reasonable approximation (especially since its dephasing effects at rate $(2T_{1})^{-1}$ can be formally included into the definition of $T_{2}$). In this limit, we can readily extremize the value of $R_{s}$ while keeping $\theta_{s}$ fixed in Eqs. (III.2), which is accomplished by fixing the ratio $\Delta_{1}^{2}/\Delta_{0}$ and then varying only $\Delta_{1}$. This procedure determines the maximum radius to be a constant, $R_{\rm max}=1/\sqrt{1/\eta+2\tau_{m}/T_{2}}<1$, that is achieved when $\Delta_{1}=\sin\theta_{s}/(\tau_{m}R_{\rm max})$. This upper bound on the radius indicates the maximum achievable purity for the set of stationary states, and thus the maximum achievable pure state preparation fidelity.

To set the control parameters to target a particular state $(y_{s},z_{s})$, we return to the $\dot{y}=\dot{z}=0$ condition for Eqs. (10) (remembering that we are dealing with the ensemble-average case where $\xi\to 0$) and solve for $\Delta_{0}$ and $\Delta_{1}$ directly to find

[TABLE]

The $\Delta_{1}$ solution is not unique, but correctly becomes unique and reduces to Eqs. (11) in the limits $T_{1},T_{2}\to\infty$, and $\eta,R_{s}\to 1$. Moreover, in the $T_{1}\to\infty$ limit the condition for the maximum radius $R_{\rm max}$ derived above corresponds precisely to unique solutions of $\Delta_{1}$ (i.e., when the square root in $\Delta_{1}$ vanishes). It thus seems reasonable to hypothesize that the condition for maximum purity corresponds to unique solutions for $\Delta_{1}$ more generally, which constrains the radius to

[TABLE]

We will revisit why this is the correct condition later when we analyze the noise-disturbance of individual trajectories in Sec. IV.2. Numerical simulations confirm that $R_{\rm max}(\theta)$ is the correct angle-dependent maximum radius that includes $T_{1}$. A plot of $R_{\rm max}(\theta)$ is shown in Fig. 5 of Sec. IV.3.

The angle-dependent $T_{1}$ corrections become significant only when $T_{1}\sim\tau_{m}$, except around the poles. That is, for sufficiently strong measurements (shorter $\tau_{m}$), the effects of $T_{1}$ can be almost safely neglected. However, for an angular range between $\delta\theta=\pm\arcsin\sqrt{2\tau_{m}/T_{1}}$ near the poles, $T_{1}$ has a significant effect: It prevents stabilization of the excited state pole, while enhancing the stability of the ground state pole (see Fig 5). Outside of this narrow angular range near each measurement pole, the limit of $\tau_{m}/T_{1}\to 0$ should be adequate to describe most experiments.

For the $T_{1}\to\infty$ limit, the required control parameters have a simple closed form in terms of only the angle $\theta$ in the $yz$-plane

[TABLE]

and generalize Eq. (11) in a straightforward way. This generalization is one of our main results, and is also shown in Fig. 1 for completeness.

III.3 Non-Markovian experimental considerations

A fundamental assumption of our analytical treatment to this point has been the Markovian nature of the feedback implicit in the Bloch Eqs. (10), which implies an infinitesimal feedback delay $T_{d}\to 0$ and thus no signal filtering. Although approximating the system time as continuous with no feedback delay yields the closed-form solutions in Eqs. (11) and (III.2) for the stabilization parameters $\Delta_{0}$ and $\Delta_{1}$, the validity of these prescriptions may not hold for more realistic experimental situations (e.g., Ref. Campagne-Ibarcq et al. (2016b)).

To address this shortcoming and more accurately model experiment, the effects of signal filtering and feedback delay within the feedback loop should be included. Signal filtering necessarily arises physically as a result of the finite bandwidth of the feedback circuitry. Feedback delay necessarily results both from the limited signal velocity in the feedback circuitry and any additional signal processing time within the feedback loop. Anticipating the need for these effects, we have already included them in the time-discrete model outlined in Eqs. (5) and (II.1).

We model the finite circuitry bandwidth as a single-pole (RC) low-pass filter with exponential time constant $T_{s}$. That is, the output filtered signal at time step $k$ is an exponential moving average $\tilde{\bar{r}}_{k}=(dt/T_{s})\sum_{j\leq k}\bar{r}_{j}\exp[-(k-j)dt/T_{s}]$ of the preceding collected raw signals $\{\bar{r}_{j}\}_{j\leq k}$. This moving average has the convenient recursive form

[TABLE]

In the limit $T_{s}\to 0$ the bandwidth becomes infinite and the signal fed back into the control loop will be identical to the measured noisy signal. As the decay constant $T_{s}$ increases, the noisy signal will become increasingly smooth and lagged compared to the raw signal. As such, we expect the quality of the stabilization to degrade with increasing $T_{s}$, since the control signal will correspond less to the instantaneous qubit state. We consider the effect of bandwidths ranging from an ideal decay constant $T_{s}\to 0$ to the measurement collapse timescale $T_{s}\to\tau_{m}$.

We model the feedback delay $T_{d}=n_{d}\,dt$ by buffering $n_{d}$ time steps $dt$ of the filtered signal before feeding that buffer back into the unitary control in Eqs. (II.1). As such, the effective Rabi frequency for the unitary control at time step $t_{k}$ is $\Delta_{0}+\Delta_{1}\,\tilde{\bar{r}}_{k-n_{d}}$. Typical values of $T_{d}$ in current laboratory setups are $\sim\!200$–$500$ns Vijay et al. (2012); Campagne-Ibarcq et al. (2016b). As noted earlier, the lag in the feedback signal due to both signal filtering and delay makes the dynamics intrinsically non-Markovian, which largely precludes better analytical estimates of the optimal stabilization parameters. Nevertheless, we numerically explore the effects below.

As shown in the left column of Fig. 2, the effect of signal filtering on the quality of the stabilized state are small. For a (representative) targeted state with angle $\theta_{T}=3\pi/10$, the stabilized angle $\theta_{s}$ drifts nearly linearly by a small amount $\pi/10$ closer to the nearest measurement pole over the large filtering range $T_{s}\in[0,\tau_{m}]$. The stabilized Bloch radius $R_{s}$ is degraded from a target value $R_{\rm max}=0.64$ by an amount $0.1$ over the same range of $T_{s}$. Notably, the radius remains essentially unaffected for $T_{s}\lesssim 0.2\tau_{m}$. Since the radius is closely related to the state purity, this robustness to small amounts of signal filtering indicates that quality state stabilization may still be achieved after a suitable recalibration of the targeted stabilization angle. Such correction of the angular drift will also compensate for most of the degradation in $R$ observed here.

As shown in the right column of Fig. 2, the effect of feedback delay on the stabilization quality is dramatic compared to the effect of signal filtering. While the (representative) target angle $\theta_{T}=3\pi/10$ drifts by approximately the same small amount $\pi/10$ toward the nearest measurement pole over the delay range $T_{d}\in[0,\tau_{m}]$, the Bloch radius $R$ becomes substantially more degraded from $R_{\rm max}=0.64$ to a mere $R_{s}=0.15$ over the same range. Such dramatic degradation of $R$ may not be mitigated by correction of the angular drift. Moreover, the Bloch radius sharply decreases linearly for even small amounts of delay, indicating that the stabilization protocol is quite sensitive to feedback delay. To be useful, this implies that the feedback delay must be a small fraction of the measurement collapse time $T_{d}/\tau_{m}\lesssim 0.2$. For typical delay times of $T_{d}\sim 200$ns, this mandates very slow measurements with $\tau_{m}\gtrsim 1\mu$s. We note that the sensitivity to delay seen here is in good qualitative agreement with a similar analysis of stabilized Rabi oscillations for a double-quantum-dot performed in Ref. Ruskov et al. (2004).

The notable exceptions to the behavior described above are the specific target angles $\theta_{T}=\pm\pi/2$ on the equator of the Bloch sphere. For these angles, stabilization is excellent when $T_{d}=T_{s}=0$, even in the non-ideal case, but rapidly becomes impossible for even small amounts of smoothing and delay. Curiously, the regions around the equator are still well-stabilized in accordance with the preceding paragraphs, but there is a sharp instability precisely at the equator. Such an instability at the equator is reminiscent of that seen in fluorescence-based qubit state stabilization protocols Wang and Wiseman (2001); Wiseman et al. (2002); Wang et al. (2005); Campagne-Ibarcq et al. (2016b); however, it is notable that our protocol does not display such an instability until signal filtering or feedback delay are added.

IV Trajectory Stabilization

In the preceding section we posed a stabilization problem for which the ensemble-average state is the only accessible quantity to the experimenter. That is, the experimenter has no knowledge of any detailed prior evolution of the state during the stabilization process, and only receives a fluctuating output state. We now consider a slightly modified posing of the problem that permits deeper investigation into the nature of the stabilization process itself.

Consider again a closed, continuous measurement feedback loop used as a storage mechanism for preparing a target qubit state on demand. However, now let the feedback loop also internally track the evolution of the qubit using the measured continuous readout. An experimenter will press a button as before, at which point the feedback loop will be terminated and the prepared qubit given to the experimenter. However, the state preparation procedure will now also report its estimate of the actual qubit state to the experimenter. The experimenter may then choose to discard the qubit and wait for a new preparation, or choose to perform an experiment immediately with the prepared state. Can the feedback parameters $\Delta_{0}$ and $\Delta_{1}$ be fixed such that a particular qubit state is reliably prepared with high probability for each single shot of the experiment?

Note that this trajectory-specific stabilization problem can dramatically differ from the ensemble-averaged version in the previous section, since the peak of the distribution of prepared states need not correspond to the mean of that distribution. The peak of the trajectory distribution corresponds to the most likely state to be prepared, which we can take to be the most relevant quantity. That is, the experimenter can readily discard the states far away from the peak while still keeping a substantial fraction of the prepared states. Thus, we seek to characterize the dominant distribution peak and the spread of the distribution around that peak.

We show that for idealized Markovian feedback, perfect stabilization of individual qubit trajectories is possible in the time-continuous limit, with an identical solution for the feedback parameters as obtained for the ensemble-averaged case. We then generalize the protocol to add realistic experimental non-idealities and show that the stabilization condition may no longer be completely satisfied for individual trajectories. Nevertheless, the optimal stabilization for individual trajectories still coincides with the stabilization of ensemble averaged states possessing maximum purity, yielding again the same solution for the feedback parameters discussed in the previous section. The lack of perfect stabilization has a substantial impact on the structure of the state trajectory distribution, however, causing broadening relative to the ideal case, as well as a shift of the most probable peak away from the ensemble-average. Finally, we discuss the effect of realistic feedback delay and signal filtering on the stabilization protocol and show that it rapidly broadens the distributional width even further unless the delay time is short compared to the measurement collapse timescale.

IV.1 Stationary states for ideal Markovian feedback

To stabilize an individual trajectory to a particular stationary point, the effect of the fluctuating noise must vanish at that point. As such, we can re-examine the Bloch Eqs. (10) in the presence of an arbitrary noise realization $\xi$ (assuming ideal limits $T_{1},\,T_{2}\to\infty$ and $\eta,\,R_{s}\to 1$ as in Sec. III.1) and demand that $\dot{y}=\dot{z}=0$ as before. Forcing the noise term to vanish for all noise realizations $\xi$ produces the additional constraints

[TABLE]

Both these constraints are automatically satisfied by the feedback parameter $\Delta_{1}=y_{s}/\tau_{m}$ already derived for the ensemble-averaged case in Eq. (11). That is, at the stationary state of the ensemble average, the backaction from the noise vanishes for each individual trajectory, so each trajectory should also be independently stabilized. This result is expected because the ensemble average solution was a pure state, indicating a certainty only obtainable if noise fluctuations identically vanish.

For time-discrete simulations, however, this ideal stability point is imperfect. The finite state jumps due to the finite time steps cause minor fluctuations around the stability point in practice, which vanish in the limit $dt\to 0$. Interestingly, near the equator with $z=0$, the width of the fluctuation distribution around the stability point with larger $dt$ narrows considerably. This improved stability at the equator is akin to balancing a rod nearly vertically by shifting its balance point laterally, which can be accomplished most efficiently when the rod is at the unstable equilibrium point at the vertical. Conversely, near the poles with $z=\pm 1$ the distribution broadens due to the destabilizing influence of the nearby pole. This converse situation is akin to balancing a rod at an extreme angle with lateral motion, which is prone to larger error since the feedback competes with the nearby stable equilibrium point. This intuition about balancing a rod also helps to explain the disastrous effect of the signal filtering and feedback delay seen at the Bloch sphere equator in Sec. III.3: If control is applied too late or in the incorrect direction for a vertically balanced rod, it will fall before the control can effectively recover. However, the stabilization of angles slightly away from the unstable equilibrium is more tolerant to delay, as the required feedback signal is already more significant.

Example trajectories are shown in Fig. 3, showing the stabilization from an initial state at angle $\theta=\pi/10$ to a target state at angle $\theta_{s}=3\pi/10$. Not only does the ensemble-average of $10^{4}$ trajectories (red) correspond to the analytical solution (blue) of the ensemble-average Markovian feedback Eq. (II.1), but most individual trajectories (gray) correctly converge to the same ensemble-average stability point within a few collapse times $\tau_{m}=0.2\mu$s. For simulation purposes, $dt=0.5$ns; note that further decreasing the time step size $dt$ will only improve the convergence to the stability point, as the analytic results were derived strictly for the $dt\to 0$ limit.

Notably, the state trajectory stability shown here behaves differently from previously considered fluorescence feedback protocols Wang and Wiseman (2001); Wang et al. (2005); Campagne-Ibarcq et al. (2016b), which were unable to stabilize the Bloch equator (even without signal filtering or feedback delay). The reason for the difference in the dispersive measurement case is that the measurement backaction symmetrically attracts the state toward both poles, so can balance the equator akin to the vertical balancing of a rod at an unstable equilibrium. In contrast, the fluorescence measurements considered in Refs. Wang and Wiseman (2001); Wiseman et al. (2002); Wang et al. (2005); Campagne-Ibarcq et al. (2016b) have fundamentally asymmetric backaction, which resulted in the equator becoming unstable.

IV.2 Minimum disturbance Markovian feedback

When nonidealities are added as in Sec. III.2, it is no longer possible to perfectly remove the backaction from the noise, even in the continuous Markovian limit $dt\to 0$. From Eq. (10) we can isolate the disturbance caused to the $y$ and $z$ coordinates per unit noise $\xi/\sqrt{\tau_{m}}$,

[TABLE]

These noise-disturbances should be minimized at a stationary point $(y_{s},z_{s})$, and depend only on the feedback parameter $\Delta_{1}$. However, they do not both vanish simultaneously for any choice of $\Delta_{1}$ if the stabilized state is not pure. That is,

[TABLE]

This discrepancy implies that a minimum fluctuation must persist for any choice of $\Delta_{1}$.

A natural resolution to this dilemma is to minimize the total noise-disturbance. That is, after defining the noise-disturbance for the entire Bloch vector $\delta\vec{s}=(0,\delta y,\delta z)$, a natural cost function is the length of this vector $C(\Delta_{1})=|\delta\vec{s}|^{2}=(\delta y)^{2}+(\delta z)^{2}$. Minimizing this cost via $C^{\prime}(\Delta_{1})=0$ yields an optimal feedback parameter that is distinct from both single-coordinate optima in Eq. (IV.2),

[TABLE]

Notably, this minimum noise-disturbance $\Delta_{1}$ precisely corresponds to the form found for the ensemble-average stability condition in Eq. (15) when the average radius $R_{s}$ is made maximal via the constraint in Eq. (III.2). We can thus understand the mysterious condition in Eq. (III.2) as physically corresponding to the requirement that the noise-disturbance be minimal for individual trajectories.

We thus completely recover the optimal results outlined for the ensemble average case in Eqs. (15) and (III.2), with limit as $T_{1}\to\infty$ given explicitly in Eq. (III.2). This result is sensible, since the ensemble-average should be most pure when each individual trajectory is optimally stabilized at its maximum achievable purity. Pragmatically, this means that there is no difference in procedure between trying to stabilize an ensemble-averaged state with maximal purity and trying to stabilize individual trajectories to have minimum measurement disturbance in the case of Markovian feedback. The two variations of our stabilization problem produce identical parameter prescriptions for $\Delta_{0}$ and $\Delta_{1}$.

Example trajectories are shown in Fig. 4, showing the same stabilization task from an initial state at angle $\theta=\pi/10$ to a target state at angle $\theta_{s}=3\pi/10$ as in Fig. 3 for comparison. Due to the nonidealities of efficiency $\eta=0.41$, energy decay time $T_{1}=60\mu$s, and environmental dephasing time $T_{2}=40\mu$s, the target maximum ensemble-averaged radius $R_{s}=0.64$ follows from Eq. (III.2), using the same collapse time $\tau_{m}=0.2\mu$s, and time step $dt=0.5$ns. As with the ideal case, the ensemble-average of $10^{4}$ trajectories (red) corresponds well to the analytic solution (blue) of the ensemble-average Markovian feedback, Eq. (II.1). However, the distribution of individual trajectories (gray) is now much broader around the ensemble-average stability point. Despite the chosen feedback parameters minimizing the total noise disturbance, there is still significant fluctuation within the stabilized region.

A steady-state trajectory histogram on the Bloch plane typically shows a single dominant lobe centered near the target stability point. Notably, however, for stabilization points near the poles at $z=\pm 1$, the distribution of trajectories may bifurcate into two distinct lobes such that the ensemble-average point corresponds to the peak of neither lobe. Example histograms of cases at $\theta_{s}=3\pi/10$ and $\theta_{s}=\pi/10$ are shown in Fig. 6 to emphasize this splitting of stabilization lobes that can occur near a pole. For the posing of the stabilization problem in this section, an experimenter will obtain likely states at the peaks of each lobe with high probability, and will recover the ensemble-average from the posing in the last section only after averaging many of these high-probability preparations.

Both the primary and secondary stability lobes for targeted angles near a pole tend to be more pure than stabilized states near the equator. However, averaging the lobes together preserves the same, nearly angle-independent, ensemble-average Bloch radius. As such, it is significantly more challenging to relate the peak of the dominant stabilized lobe (i.e., the most likely prepared state) to the two control parameters $\Delta_{0}$ and $\Delta_{1}$ analytically. Nevertheless, the positions of the peaks for the dominant lobes may be readily mapped out numerically by varying the ostensible target angles $\theta_{s}$ according to Eqs. (22) and (III.2) and computing the obtained $\theta_{P}$ and $R_{P}$ for the dominant histogram peak corresponding to those parameters. We plot the result of this numerical procedure in Fig. 5, showing the purification of the most likely states near the measurement poles.

IV.3 Non-Markovian experimental considerations

As with the ensemble average posing in Sec. III.3, adding single-pole low-pass signal filtering and feedback delay dramatically alters the quality of the stabilization protocol. The effect can be seen most readily in the trajectory histograms, of which we show an example in Fig. 7 that can be directly compared with the top plot in Fig. 6. The single lobe stabilization region spreads out due to both signal filtering and feedback delay into a larger region, consistent with the expectations from Fig. 2. As such, while most prepared states will actually have relatively high purity, the ensemble average state becomes increasingly mixed. Moreover, the angular stabilization is compromised as the angular uncertainty in the broadened region becomes large. While the most probable state remains near the targeted angle, with high purity, it becomes comparable in probability to most other prepared states, and thus dubiously stabilized.

In order to achieve the best quality stabilization, an experimenter should thus minimize the feedback delay $T_{d}$ relative to the measurement collapse time $\tau_{m}$. As estimated in Sec. III.3, for realistic delays of $T_{d}\sim\!200$ ns, one should use slow measurements with $\tau_{m}\gtrsim 1$ $\mu$s. Compared to this timescale, typical signal filtering timescales $T_{s}$ should also be relatively inconsequential.

V Conclusions

We investigated the linear feedback stabilization of a transmon qubit being continuously monitored by dispersive coupling to a microwave field. Using only a time-varying Rabi drive consisting of a constant frequency modulated by a direct linear feedback of the collected and filtered microwave signal, we showed that it is possible to prepare target qubit states on demand for any angle in the Bloch sphere. We detailed both a time-discrete and a time-continuous model for the measurement process, derived analytical solutions for the stabilization conditions, and checked the results through numerical simulations. Notably, we showed that the required feedback parameters are identical for two distinct variations of the stabilization problem: (1) on-demand preparation of a stabilized ensemble-average state, and (2) on-demand preparation of a single state trajectory with minimized total measurement disturbance. In the first situation, an experimenter does not have any information about the exact preparation shot-to-shot and so can only access the ensemble-averaged state. In the second, the experimenter also has information about the exact state that was prepared. Interestingly, in the latter the most likely states that are prepared tend to have higher purity, yet still recover the same ensemble average state as the former.

In the ideal case we showed that any target pure qubit state may be stabilized and derived the two unique feedback parameters in Eq. (11) that achieve that stabilization. In the nonideal case, with experimental nonidealities of energy relaxation, environmental dephasing, and measurement inefficiency, we showed how to target ensemble average states with maximum purity by setting the feedback parameters in Eqs. (22) and (15), which have the approximate closed form in Eq. (III.2) when energy relaxation can be neglected. Expressed in Bloch coordinates, most qubit state angles may still be stabilized in the nonideal case, with a nearly angle-independent maximum radius given by Eq. (III.2). Finally, we included the experimental nonidealities of single-pole low-pass signal filtering and feedback delay, and investigated their effect numerically. These additional effects physically arise from the feedback circuitry, and degrade the quality of the achievable state stabilization by widening the distribution of prepared state trajectories and thus shrinking the ensemble averaged state purity. Notably, feedback delay that is appreciable compared to the collapse timescale has the most damaging effect on the stabilization protocol.

This stabilization protocol is within reach of modern superconducting transmon experiments involving continuous microwave measurements. The analysis presented here should be suitable for accurately modeling a feedback loop implemented with a simple controller programmed into a field-programmable gate array (FPGA) as used in Refs. Campagne-Ibarcq et al. (2013, 2016a, 2016b), or even a feedback loop consisting of a simple analog circuit as used in Refs. Vijay et al. (2012); de Lange et al. (2014). Such a feedback stabilization experiment would prompt further investigation into the use of continuous measurement feedback for practical quantum information tasks in future work.

Acknowledgements.

We thank Alexander Korotkov, Howard Wiseman, Rusko Ruskov, Joshua Combes, Shay Hacohen-Gourgy, Leigh Martin, Philippe Lewalle, and Irfan Siddiqi for helpful discussions. This work was supported by U.S. Army Research Office Grants No. W911NF-15-1-0496 and No. W911NF-13-1-0402, by National Science Foundation Grant No. DMR-1506081 including an REU extension, and by the Development and Promotion of Science and Technology Talents Project Thailand. We also acknowledge partial support by the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development and Innovation.

Appendix

In this appendix, we briefly include a pure state formulation of the continuous time limit for completeness. This formulation simply reproduces the conclusions already obtained from the Bloch equations in the main text. In the ideal case without experimental nonidealities, Eqs. (3) and (4) are sufficient to describe the pure state vector update equation

[TABLE]

Performing an expansion to linear order in $dt$ of the pure state update Eq. (23) while applying the Itô rule $dW^{2}=dt$ then produces the Itô stochastic Schrödinger equation

[TABLE]

Though this state vector representation has lower dimensionality than the density operator SME in Eq. (II.2), it reproduces the same Bloch coordinate equations in Eq. (10), albeit with the added constraint $x^{2}+y^{2}+z^{2}=1$.