# Stochastic modeling of superconducting qudits in the dispersive regime

Kangdi Yu<EMAIL_ADDRESS>Murat C. Sarihan  Jin Ho Kang  Madeline Taylor
Cody S. Fan  Ananyo Banerjee  Department of Electrical and Computer
Engineering, University of California, Los Angeles, California 90095, USA
Jonathan L. DuBois  Yaniv J. Rosen  Lawrence Livermore National Laboratory,
Livermore, California 94550, USA  Chee Wei Wong<EMAIL_ADDRESS>Department of Electrical and Computer Engineering, University of California,
Los Angeles, California 90095, USA

###### Abstract

The field of superconducting quantum computing, based on Josephson junctions,
has recently seen remarkable strides in scaling the number of logical qubits.
In particular, the fidelities of one- and two-qubit gates have reached the
breakeven point with the novel error mitigation and correction methods.
Parallel to these advances is the effort to expand the Hilbert space within a
single junction or device by employing high-dimensional qubits, otherwise
known as qudits. Research has demonstrated the possibility of driving higher-
order transitions in a transmon or designing innovative multimode
superconducting circuits, termed multimons. These advances can significantly
expand the computational basis while simplifying the interconnects in a large-
scale quantum processor. In this work we extend the measurement theory of a
conventional superconducting qubit to that of a qudit, focusing on modeling
the dispersive quadrature measurement in an open quantum system. Under the
Markov assumption, the qudit Lindblad and stochastic master equations are
formulated and analyzed; in addition, both the ensemble-averaged and the
quantum-jump approach of decoherence analysis are detailed with analytical and
numerical comparisons. We verify our stochastic model with a series of
experimental results on a transmon-type qutrit, verifying the validity of our
high-dimensional formalism.

††preprint: APS/123-QED

##  I Introduction

Superconducting quantum computation based on single- and two-qubit gates [1,
2, 3] has recently demonstrated milestone successes in comparison to classical
computation [4, 5]. While the error correction scheme can overcome the noise
in the qubit-qubit interaction and correct the undesired decoherence due to
coupling to the environment, the required hardware resources, chip footprint,
and peripheral connections will soon face scalability issues. To create a more
hardware-efficient quantum process with fewer quantum gates and interconnects,
an emerging effort has been placed into developing high-dimensional
superconducting quantum computing. In particular, the hope is that both the
size of the quantum algorithm and the error rate induced in a long gate
sequence can be reduced by using the naturally available higher energy levels
of the localized artificial atoms, also termed qudits.

Experiments involving the usage of higher energy levels already in the
transmon qubits have been examined [6, 7]. Furthermore, a more complex
Josephson-junction-based network with multiple nonlinear modes coupled
longitudinally [8, 9] can also be developed. Despite the recent experimental
efforts on higher dimensional quantum computation, the theory of dispersive
measurement used in inferring a qudit state under sources of decoherence still
mostly relies on a heuristic extension from the qubit formalism. In this work,
we adopt both Lindblad master equations and the method of quantum trajectories
to analyze the state evolution of a qudit when measured dispersively.

In circuit or cavity quantum electrodynamics (QED), continuous measurement and
other unwanted coupling to the environment can be modeled by master equations
of the Lindblad form under the Born-Markov assumption. Analysis of the qubit-
resonator coupling in the dispersive regime has relied on this formalism,
aided by the good agreement between measurement and the predicted measurement-
induced dephasing and number splitting [10, 11]. Our work first generalizes
the measurement of the observable $\hat{\sigma}_{z}$ in the qubit case to a
generic measurement in the longitudinal direction of the qudit, thus extending
the notion of measurement-induced dephasing. Due to the longitudinal coupling
appearing in the dispersive regime, we show analytically that each energy
eigenstate of the qudit is entangled with a coherent state of the resonator;
consequently, measurement of the quadrature fields of the resonator can be
used to infer the qudit state.

However, since the qudit Lindblad master equation only describes the ensemble-
averaged time evolution of the state, it does not make use of the information
leaking out of the system, part of which is, of course, what one measures in
real experiments. To capture the update of an observer’s knowledge during a
continuous measurement, in this work we subsequently apply the quantum
trajectory theory [12, 13, 14, 15] in which our knowledge of the quantum state
can be modified based on the measurement record. Given that the measurement
outcomes are stochastic as required by quantum mechanics, the entire
measurement record forms a random process that can be modeled as a diffusive
process when the measurement is continuous and weak. In particular, under the
diffusive limit, we derive an effective heterodyne stochastic master equation
(SME) of the qudit subspace (in the Itô sense) which is then compared with the
unconditioned master equation after taking the ensemble average. Recently, the
stochastic nature of qudit dispersive measurement has been theoretically
studied in the absence of other decoherence processes such as qudit
$T_{1}$-decay and dephasing, and specific examples such as the $N$-level
“clock” system were examined [16]. Our approach relaxes the assumptions to
include other experimentally observed decoherence processes in a transmon
qudit and makes an intuitive link between the stochastic trajectories solved
from SME and the averaged dynamics given by the Lindblad master equation.

In this work, we first review the dispersive coupling of a high-dimensional
system to a linear resonator in Section II. In Section III, the
(unconditioned) Lindblad master equation of the combined system (i.e., a qudit
plus a resonator) is solved analytically using two methods – the positive
$P$-representation and the qudit-state-dependent displacement operator.
Consequently, an effective master equation for the qudit with the resonator
degrees of freedom traced out can be fortunately formulated as a Markovian
system when the measurement-induced frequency shifts are negligible. Then, in
Section IV, we turn to the conditioned qudit state and derive an effective
qudit SME in the diffusive limit, with the measurement outcomes described by
two stochastic differential equations. Finally, in Section V, the simulated
quantum trajectories are compared with our experiments on a transmon qutrit
(i.e., a three-level system) coupled to a 3D cavity, demonstrating a
remarkable match between our formalism and experiments.

##  II cQED for Dispersive Measurement of a Qudit

For a general analysis, we consider a weakly anharmonic oscillator coupled to
a linear resonator for readout and control. For future reference, we denote
the anharmonic oscillator, also referred to as the system or qudit, by
$\mathcal{S}$ and the resonator by $\mathcal{R}$. Hence, the combined system
lives in the Hilbert space
$\mathscr{H}_{\mathcal{SR}}=\mathscr{H}_{\mathcal{S}}\otimes\mathscr{H}_{\mathcal{R}}$.
The anharmonic oscillator model can be used to describe a transmon-type
superconducting qudit while the resonator, modeled by a harmonic oscillator,
can be realized as a 3D cavity or a planar transmission-line resonator [17,
18, 19]. On the one hand, the size of a typical 3D cavity, and thus the
associated mode volume, is much larger than that of a superconducting qudit,
allowing one to apply the dipole approximation to describe the qudit-resonator
interaction. On the other hand, even though a planar resonator has a more
confined mode profile and is coupled to the qubit locally (e.g., capacitive
coupling), one can still derive a similar dipole interaction at the circuit
level. Consequently, for both the 3D and 2D resonators, the Hamiltonian of the
combined system under the rotating-frame approximation (RWA) is given by [20]

$\hat{H}/\hbar=\omega_{\text{q}}\hat{a}_{\text{q}}^{\dagger}\hat{a}_{\text{q}}+\frac{\alpha_{\text{q}}}{2}\hat{a}_{\text{q}}^{\dagger}\hat{a}_{\text{q}}^{\dagger}\hat{a}_{\text{q}}\hat{a}_{\text{q}}\\\
+\omega_{\text{r}}\\!\left(\hat{a}_{\text{r}}^{\dagger}\hat{a}_{\text{r}}+\frac{1}{2}\right)-\Big{(}g\hat{a}_{\text{r}}\hat{a}_{\text{q}}^{\dagger}+g^{*}\hat{a}_{\text{r}}^{\dagger}\hat{a}_{\text{q}}\Big{)},$
(1)

where $\omega_{\text{q}}$ is the qubit frequency (i.e., the transition
frequency between the ground and first excited states), $\alpha_{\text{q}}$
the anharmonicity of the qudit, $\omega_{\text{r}}$ the resonator frequency,
and $g$ the qudit-resonator coupling coefficient. We also define
$\Delta_{\text{qr}}=\omega_{q}-\omega_{r}$ as the detuning between the qudit
and the resonator.

In Eq.(1), we have assumed the fourth-order expansion of the transmon
Hamiltonian, which ignores the fact that a transmon can only support a finite
number of bound energy eigenstates. More generally, we can replace the fourth-
order expansion with the Hamiltonian of a general qudit such that

$\hat{H}/\hbar=\sum_{j=0}^{D-1}\omega_{j}\ket{j}\\!\bra{j}+\omega_{\text{r}}\\!\left(\hat{a}_{\text{r}}^{\dagger}\hat{a}_{\text{r}}+\frac{1}{2}\right)\\\
-\sum_{j,k=0}^{D-1}\Big{(}g_{jk}\ket{j}\\!\bra{k}\hat{a}_{\text{q}}^{\dagger}+g_{jk}\ket{k}\\!\bra{j}\hat{a}_{\text{q}}\Big{)}$
(2)

with $\omega_{j}$ and $\ket{j}$ representing the energy and state vector of
the $j$th energy level of the qudit. The weakly anharmonic model corresponds
to the case where

$g_{jk}\approx\begin{cases}\sqrt{j+1}\,g&\text{if }\ \ j-k=1,\\\
0&\text{otherwise},\end{cases}$ (3)

for all $j,k=0,...,D-1$ ($g$ is the same coupling coefficient defined in
Eq.(1)).

To perform a quantum non-demolition (QND) measurement (justified in the
following text), we set the resonator frequency to be detuned from (and
usually higher than) the qubit frequency. In the dispersive regime where
$|g|\ll|\Delta_{\text{qr}}|$, Eq.(1) can be approximated as

$\hat{H}^{\text{disp}}/\hbar=\tilde{\omega}_{\text{q}}\hat{a}_{\text{q}}^{\dagger}\hat{a}_{\text{q}}+\frac{\alpha_{\text{q}}}{2}\hat{a}_{\text{q}}^{\dagger}\hat{a}_{\text{q}}^{\dagger}\hat{a}_{\text{q}}\hat{a}_{\text{q}}\\\
+\omega_{\text{r}}\\!\left(\hat{a}_{\text{r}}^{\dagger}\hat{a}_{\text{r}}+\frac{1}{2}\right)+\chi_{\text{qr}}\hat{a}_{\text{q}}^{\dagger}\hat{a}_{\text{q}}\hat{a}_{\text{r}}^{\dagger}\hat{a}_{\text{r}},$
(4)

where $\tilde{\omega}_{\text{q}}=\omega_{\text{q}}+|g|^{2}/\Delta_{\text{qr}}$
is the Lamb-shifted qubit frequency and

$\chi_{\text{qr}}=\frac{2\alpha_{\text{q}}|g|^{2}}{\Delta_{\text{qr}}(\Delta_{\text{qr}}+\alpha_{\text{q}})}$
(5)

is the dispersive shift, also known as the (fourth-order) cross-Kerr
coefficient [20, 21]. By lumping the last term in Eq.(4) into the resonator
Hamiltonian, one observes a qudit-state-dependent shift in the resonator
frequency. In particular, if the qudit is in the energy eigenstate $\ket{j}$,
the resonator will experience a dispersive shift $j\chi_{\text{qr}}$. By
determining this frequency shift via a resonator transmission or reflection
measurement, one should be able to infer the qudit state.

Moreover, if Eq.(2) is adopted instead, one finds

$\hat{H}^{\text{disp}}/\hbar=\sum_{j=0}^{D-1}(\omega_{j}+\Lambda_{j})\ket{j}\\!\bra{j}\\\
+\omega_{\text{r}}\\!\left(\hat{a}_{\text{r}}^{\dagger}\hat{a}_{\text{r}}+\frac{1}{2}\right)+\sum_{j=0}^{D-1}\chi_{j}\hat{a}^{\dagger}\hat{a}\ket{j}\\!\bra{j}.$
(6)

to the second order in $|g_{jk}/(\omega_{j}-\omega_{k}-\omega_{\text{r}})|$ in
the dispersive regime. The Lamb shift $\Lambda_{j}$ and the dispersive shift
$\chi_{j}$ of the $j$th energy level of the qudit are given, respectively, by
[22]

$\Lambda_{j}=\sum_{k=0}^{D-1}\chi_{jk}=\sum_{k=0}^{D-1}\frac{|g_{jk}|^{2}}{\omega_{j}-\omega_{k}-\omega_{\text{r}}},$
(7) $\chi_{j}=\sum_{k=0}^{D-1}(\chi_{jk}-\chi_{kj}),$ (8)

where

$\chi_{jk}=\frac{|g_{jk}|^{2}}{\omega_{j}-\omega_{k}-\omega_{\text{r}}}.$ (9)

The main goal of this work is to quantify the qudit dispersive measurement in
terms of the rate at which the information leaks out from the resonator. In
addition, we need to answer to what extent a weak and continuous dispersive
measurement is QND. The following analysis does not rely on the particular
form of the dispersive shift nor on the relationships among the dispersive
shifts of different energy levels; for simplicity, we will mainly use Eq.(4).
Nevertheless, it should be mentioned that Eq.(4) and (6) are valid only when
the resonator photon number is low [22]. In particular, the Schrieffer–Wolff
transformation used in deriving Eq.(4) and (6) assumes that the strength of
the qudit-resonator interaction is much smaller than $\Delta_{\text{qr}}$. Due
to the creation and annihilation operators that appear in the qudit-resonator
coupling term, there is a $\sqrt{n_{\text{r}}}$-scaling of the interaction
strength, where $n_{\text{r}}$ is the resonator photon number. Hence, we
constrain to a low readout power in both the derivation and experiment.

##  III Unconditioned Master Equation

We first consider the unconditioned master equation where information leaking
out from the combined system is averaged unconditionally. Later, we derive a
stochastic differential equation in which the density operator at time $t$ is
conditioned on the heterodyne measurement at time $s<t$. Since only the
ensemble-averaged state is examined in the unconditioned master equation, the
solution to the differential equation is deterministic given an initial
condition.

The master equation of a qubit coupled to a resonator dispersively has been
studied extensively [10, 23]. It has been shown that each energy eigenstate
(i.e., the $z$-basis vectors) of the qubit is entangled with a coherent state
of the resonator; in addition, a qubit in a superposition state subject to a
continuous readout pulse will dephase without changing the population in the
$z$-basis. Here we generalize the solution to an arbitrary qudit in the
dispersive-coupling regime. To avoid writing down too many equations when $D$
is large, we will explicitly show the derivation for a qutrit measured
dispersively while experiencing the $T_{1}$-decay and pure dephasing;
nevertheless, the approaches used for solving the master equation can be
easily extended to higher dimensional systems by adding terms of similar
forms.

To set up the problem, we let the combined system be a qutrit (labeled as
$\mathcal{S}$) coupled to a resonator ($\mathcal{R}$) dispersively. Note that
the environment is not a part of the composite system, i.e., we have already
traced out the environment to write down a master equation. The (reduced)
state of the combined system, denoted by $\hat{\rho}_{\mathcal{SR}}$, lives in
the Hilbert space $\mathscr{H}_{\mathcal{SR}}$. For clarity, we denote the
energy levels of the qutrit by $\ket{g}$, $\ket{e}$, and $\ket{f}$ (ordered
with increasing energy) to differentiate from the Fock states of the
resonator. We study the time evolution of the composite state under the usual
Born and Markov approximations [24, 25] in which the state transition maps
form a quantum dynamical semigroup and are described by a Lindblad master
equation.

###  III.1 Master Equation for the Composite System in the Laboratory Frame

Suppose the resonator is coupled to the environment with a total decay rate of
$\kappa$. For a 3D microwave cavity with two ports, $\kappa$ is the sum of the
input decay rate $\kappa_{\text{in}}$, output decay rate
$\kappa_{\text{out}}$, and the internal decay rates $\kappa_{\text{int}}$ due
to material losses. If the resonator is configured in the reflection mode
(i.e., one port), then
$\kappa_{\text{in}}=\kappa_{\text{out}}\doteq\kappa_{\text{ext}}$ and the
total decay rate is $\kappa=\kappa_{\text{int}}+\kappa_{\text{ext}}$. In
reality, the resonator supports many modes; here we focus on only one mode
(usually the fundamental mode) of the resonator with frequency
$\omega_{\text{r}}$ which captures the dispersive measurement succinctly. The
resonator-environment interaction is modeled as a harmonic oscillator coupled
to a continuum of bath oscillators. At superconducting temperature, we assume
that the bath is in the vacuum state (i.e., the mean photon number at the
resonator frequency is $\bar{N}(\omega_{\text{r}})=0$) so that the usual terms
in the quantum optical master equation [26], i.e.,

$\displaystyle\kappa\big{[}\bar{N}(\omega_{\text{r}})+1\big{]}\mathcalboondox{D}\big{[}\hat{a}\big{]}\hat{\rho}_{\mathcal{R}}(t)+\kappa\bar{N}(\omega_{\text{r}})\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{a}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{R}}(t)$

reduces to $\mathcalboondox{D}\big{[}\hat{a}\big{]}\hat{\rho}_{\mathcal{R}}$,
where
$\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}$
is the dissipation superoperator associated with the Lindblad operator
$\hat{L}$ defined via the action

$\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}=\hat{L}\hat{\rho}\hat{L}^{\dagger}-\frac{1}{2}\hat{L}^{\dagger}\hat{L}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{L}^{\dagger}\hat{L}$
(10)

on any density operator $\hat{\rho}$ [27, 28, 13, 26].

For the qutrit, we study both spontaneous decay and pure dephasing. Without
imposing any selection rule, we assume the qutrit can decay from $\ket{f}$ to
$\ket{e}$, from $\ket{f}$ to $\ket{g}$, and from $\ket{e}$ to $\ket{g}$ with
decay rates $\gamma_{1,ef}$, $\gamma_{1,gf}$, and $\gamma_{1,ge}$,
respectively. We also include the pairwise pure dephasing with rates
$\gamma_{\phi,ge},\gamma_{\phi,gf},$ and $\gamma_{\phi,ef}$ to study the
coherence time of superposition states. One can show that the three pairwise
dephasing terms are equivalent to a single dephasing term with three energy
levels included; that is, given any set
$\\{\gamma_{\phi,g},\gamma_{\phi,e},\gamma_{\phi,f}\\}$, there exists
$\\{\gamma_{\phi,ge},\gamma_{\phi,gf},\gamma_{\phi,ef}\\}$ such that

$\mathcalboondox{D}\left[\sum_{a\in\\{g,e,f\\}}\sqrt{\gamma_{\phi,a}}\ket{a}\\!\bra{a}\right]\\\
=\sum_{a\in\\{g,e,f\\}}\sum_{b>a}\mathcalboondox{D}\left[\sqrt{\frac{\gamma_{\phi,ab}}{2}}\left(\ket{a}\\!\bra{a}-\ket{b}\\!\bra{b}\right)\right].$
(13)

Therefore, we stick with the three pairwise dephasing rates without loss of
generality. For a general qudit with $D$ levels, we will need $D(D-1)/2$
pairwise dephasing terms; Eq.(13) can also be extended to a $D$-level system
easily.

By including the decoherence channels mentioned above, we can write down the
Lindblad master equation of the composite system

$\displaystyle\dot{\hat{\rho}}_{\mathcal{SR}}(t)$
$\displaystyle=-\frac{\mathrm{i}}{\hbar}\\!\left[\hat{H}_{\text{eff}}(t),\hat{\rho}_{\mathcal{SR}}(t)\right]+\kappa\mathcalboondox{D}\big{[}\hat{a}\big{]}\hat{\rho}_{\mathcal{SR}}(t)+\gamma_{1,ge}\mathcalboondox{D}\big{[}\hat{\sigma}_{ge}\big{]}\hat{\rho}_{\mathcal{SR}}(t)+\gamma_{1,gf}\mathcalboondox{D}\big{[}\hat{\sigma}_{gf}\big{]}\hat{\rho}_{\mathcal{SR}}(t)+\gamma_{1,ef}\mathcalboondox{D}\big{[}\hat{\sigma}_{ef}\big{]}\hat{\rho}_{\mathcal{SR}}(t)$
$\displaystyle\ \ \ \ \ \ \
+\frac{\gamma_{\phi,ge}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,ge}\big{]}\hat{\rho}_{\mathcal{SR}}(t)+\frac{\gamma_{\phi,gf}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,gf}\big{]}\hat{\rho}_{\mathcal{SR}}(t)+\frac{\gamma_{\phi,ef}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,ef}\big{]}\hat{\rho}_{\mathcal{SR}}(t)$
$\displaystyle\approx-\frac{\mathrm{i}}{\hbar}\\!\left[\hat{H}_{\text{eff}}(t),\hat{\rho}_{\mathcal{SR}}(t)\right]+\left(\kappa\mathcalboondox{D}\big{[}\hat{a}\big{]}+\frac{\gamma_{2,ge}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,ge}\big{]}+\frac{\gamma_{2,gf}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,gf}\big{]}+\frac{\gamma_{2,ef}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,ef}\big{]}\right)\hat{\rho}_{\mathcal{SR}}(t),$
(14)

where we adopt the notations

$\displaystyle\hat{\sigma}_{z,ab}=\ket{a}\\!\bra{a}-\ket{b}\\!\bra{b}\ \
\text{ and }\ \ \hat{\sigma}_{ab}=\ket{a}\\!\bra{b}$ (15)

for $a\neq b$. In addition, by assuming that $T_{1,ab}=1/\gamma_{1,ab}$ is
much longer than the other decoherence timescales, we have removed the qutrit
decay terms and have lumped the extra dephasing rates $\gamma_{1,ab}/2$ with
the pure dephasing rates $\gamma_{\phi,ab}$ to define
$\gamma_{2,ab}=\gamma_{\phi,ab}+\gamma_{1,ab}/2$ for
$(a,b)\in\\{(g,e),(g,f),(e,f)\\}$. In the Appendices, two approaches to
solving the master equation are discussed: The first method, which requires
less algebra, is the method of positive $P$-representation, typically used
when a harmonic oscillator is driven into a coherent state [10]. The
long-$T_{1}$ assumption is used for the derivation to have a closed-form
solution. In contrast, the second method, introduced in [12], relies on the
resonator displacement operator and allows one to solve Eq.(14) with
$\gamma_{1,ab}$ included.

The effective Hamiltonian of a qutrit coupled with a resonator in the
dispersive regime subject to a readout probe (under the RWA) is given by

$\displaystyle\hat{H}^{\text{disp}}_{\mathcal{SR}}(t)/\hbar$
$\displaystyle=\tilde{\omega}_{\text{q}}\ket{e}\\!\bra{e}+(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\ket{f}\\!\bra{f}+\omega_{\text{r}}\hat{a}^{\dagger}\hat{a}$
$\displaystyle\ \ \ \
+\chi_{\text{qr}}(\ket{e}\\!\bra{e}+2\ket{f}\\!\bra{f})\hat{a}^{\dagger}\hat{a}$
$\displaystyle\ \ \ \
-\left(\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}e^{-\mathrm{i}\omega_{\text{d}}t}\,\hat{a}^{\dagger}+\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}^{*}e^{\mathrm{i}\omega_{\text{d}}t}\,\hat{a}\right),$
(16)

where we have set the zero-energy reference to be the ground-state energy of
the dressed system and used
$\tilde{\omega}_{\text{q}}=\tilde{\omega}_{e}-\tilde{\omega}_{g}$ to denote
the qubit frequency with the Lamb shift included. To address the second-
excited state $\ket{f}$, we also introduce the anharmonicity
$\alpha_{\text{q}}=(\tilde{\omega}_{f}-\tilde{\omega}_{e})-2\tilde{\omega}_{\text{q}}$,
which is negative for a transmon. The last term in Eq.(16) represents the
readout signal sent to the resonator with frequency $\omega_{\text{d}}$ and
amplitude (in terms of the square root of the average photon flux)
$\bar{a}_{\text{in}}$. It is assumed that the readout signal can be modeled by
a classical signal [17, 13] (i.e., the stiff-pump limit) and the imaginary
frequency added to the resonator due to the decay is negligible since
$\kappa\ll\omega_{\text{r}}$. Moreover, for a transmon-type qudit, we use the
fact that the dispersive shift is a linear function of the number of
excitations in the qudit, i.e., the cavity frequency shifts by
$\chi_{\text{qr}}$ when exited from $\ket{g}$ to $\ket{e}$ and shifts by
$2\chi_{\text{qr}}$ when exited from $\ket{g}$ to $\ket{f}$. We emphasize that
the two methods used to solve the master equation are still valid even if the
dispersive shift scales nonlinearly.

The subsequent calculation can be simplified if we move the cavity part of the
Hamiltonian to a frame that rotates at the drive frequency
$\omega_{\text{d}}$. (Note that the qutrit Hamiltonian stays the same.) Then,
the time-varying drive $\varepsilon_{\text{d}}(t)$ reduces to a complex scalar
$\epsilon=\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}$ and the Hamiltonian in
this rotating frame

$\displaystyle\hat{H}_{\mathcal{SR},\text{rot}}^{\text{disp}}/\hbar$
$\displaystyle=\tilde{\omega}_{\text{q}}\ket{e}\\!\bra{e}+(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\ket{f}\\!\bra{f}+\Delta_{\text{rd}}\hat{a}^{\dagger}\hat{a}$
$\displaystyle\ \ \ \
+\chi_{\text{qr}}(\ket{e}\\!\bra{e}+2\ket{f}\\!\bra{f})\hat{a}^{\dagger}\hat{a}-\left(\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a}\right),$
(17)

with $\Delta_{\text{rd}}=\omega_{\text{r}}-\omega_{\text{d}}$, is now time-
independent. Then, the master equation of the composite system in the rotating
frame is obtained by making the substitution
$\hat{H}_{\text{eff}}=\hat{H}_{\mathcal{SR},\text{rot}}^{\text{disp}}$ in
Eq.(14).

Figure 1: Schematic of qutrit readout and solution of the composite-system
master equation. a. The input-output perspective of the transmission-mode
dispersive measurement. The readout signal $\hat{a}_{\text{in}}$ entering the
cavity from the left port (i.e., port 1) is approximated by a classical drive
with complex amplitude $\bar{a}_{\text{in}}$ while the transmitted signal at
the right port (i.e., port 2) described by the traveling-wave annihilation
operator $\hat{a}_{\text{out},2}$ in order to capture the quadrature
uncertainty. b. The transient complex amplitude of the three coherent states
$\ket{\alpha_{a}}$ of the resonator associated with the $\ket{a}$ for
$a=g,e,f$. The steady state of each coherent state amplitude lies on a circle
going through the origin of the phase plane. Inset: The build-up of mean
photon number of $\ket{\alpha}_{a}$ as a function of time. c. Distance between
two coherent state amplitudes as a function of the readout frequency. To
illustrate a more general trend, we also include the fourth energy level
$\ket{h}$ of the transmon. d. Same as c, but plotted with $\chi_{\text{qr}}$
smaller, equal, or larger than a fixed $\kappa$. e/f. Time evolution of the
composite state as solved from the composite-system master equation in the
long-$T_{1}$ limit.

###  III.2 Analytical Solution to the Master Equation of the Combined System

To solve Eq.(14), we first express $\hat{\rho}_{\mathcal{SR}}$ as

$\hat{\rho}_{\mathcal{SR}}(t)=\sum_{a,b\in\\{g,e,f\\}}\hat{\rho}_{ab}(t)\ket{a}\\!\bra{b},$
(18)

where

$\hat{\rho}_{ab}(t)=\bra{a}\hat{\rho}_{\mathcal{SR}}(t)\ket{b}$ (19)

for $a,b\in\\{g,e,f\\}$ are operators acting on the Fock space of the
resonator. Now, Eq.(14) can be rewritten as a set of nine coupled operator
equations in terms of $\hat{\rho}_{ab}$. As shown in Appendix A.1, by
expressing $\hat{\rho}_{ab}$ in terms of the positive $P$-representation [29,
13]

$\hat{\rho}_{ab}(t)=\int\mathrm{d}^{2}\alpha\int\mathrm{d}^{2}\beta\frac{\ket{\alpha}\\!\bra{\beta^{*}}}{\innerproduct{\beta^{*}}{\alpha}}P_{ab}(\alpha,\beta,t),$
(20)

one can reduce the nine operator differential equations into nine scalar
equations of $P_{ab}$. The technique of the positive $P$-representation has
already been applied to the qubit case to reveal the measurement-induced
dephasing on the qubit; Appendix A demonstrates the effectiveness of this
method in solving the qutrit-resonator master equation. In fact, one can apply
it to a general qudit system (in the long-$T_{1}$ limit).

For the qutrit case, the time evolution of the combined state is shown to be

$\displaystyle\hat{\rho}_{\mathcal{SR}}(t)=\sum_{a\in\\{g,e,f\\}}p_{a}(0)\ket{a}\\!\bra{a}\otimes\ket{\alpha_{a}(t)}\\!\bra{\alpha_{a}(t)}$
$\displaystyle+\sum_{a,b\in\\{g,e,f\\}}\frac{c_{ab}(t)}{\bra{\alpha_{b}(t)}\ket{\alpha_{a}(t)}}\ket{a}\\!\bra{b}\otimes\ket{\alpha_{a}(t)}\\!\bra{\alpha_{b}(t)},$
(21)

where $p_{g,e,f}(0)$ are the initial populations in $\ket{g}$, $\ket{e}$, and
$\ket{f}$, respectively and $\ket{\alpha_{a}(t)}$ represent the coherent
states of the resonator. In addition,

$\displaystyle\dot{\alpha}_{g}$
$\displaystyle=-\mathrm{i}(\Delta_{\text{rd}}-\mathrm{i}\kappa/2)\alpha_{g}+\mathrm{i}\epsilon,$
(22) $\displaystyle\dot{\alpha}_{e}$
$\displaystyle=-\mathrm{i}(\Delta_{\text{rd}}+\chi_{\text{qr}}-\mathrm{i}\kappa/2)\alpha_{e}+\mathrm{i}\epsilon,$
(23) $\displaystyle\dot{\alpha}_{f}$
$\displaystyle=-\mathrm{i}(\Delta_{\text{rd}}+2\chi_{\text{qr}}-\mathrm{i}\kappa/2)\alpha_{f}+\mathrm{i}\epsilon.$
(24)

determine the time evolution of the coherent states entangled with the qutrit
states whereas

$\displaystyle\dot{c}_{ge}$
$\displaystyle=\mathrm{i}(\tilde{\omega}_{\text{q}}+\mathrm{i}\gamma_{2,ge})c_{ge}+\mathrm{i}\chi_{\text{qr}}\alpha_{g}\alpha_{e}^{*}c_{ge},$
(25) $\displaystyle\dot{c}_{gf}$
$\displaystyle=\mathrm{i}(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}}+\mathrm{i}\gamma_{2,gf})c_{gf}+\mathrm{i}2\chi_{\text{qr}}\alpha_{g}\alpha_{f}^{*}c_{gf},$
(26) $\displaystyle\dot{c}_{ef}$
$\displaystyle=\mathrm{i}(\tilde{\omega}_{\text{q}}+\alpha_{\text{q}}+\mathrm{i}\gamma_{2,ef})c_{ef}+\mathrm{i}\chi_{\text{qr}}\alpha_{e}\alpha_{f}^{*}c_{ef},$
(27)

with $c_{ab}=c_{ba}^{*}$ for $a\neq b$, govern the oscillation and decay of
the off-diagonal terms of the density operator. Since we have ignored
$\gamma_{1,ab}$, we observe that the populations of the qutrit do not change
over time, a critical feature of the QND measurement. However, the coherence
terms will decay to zero with additional rates proportional to
$\chi_{\text{qr}}$. It should be noted that $\alpha_{a}\alpha_{b}^{*}$ is
complex in general; hence, the last term in Eq.(25), (26) or (27) contains
both a decay and a frequency shift. An example of the time evolution of
$\hat{\rho}_{\mathcal{SR}}$ is plotted in Figure 1(e) and (f).

Given the general solution, of particular interest are the steady-state
amplitudes of the cavity coherent states

$\displaystyle\alpha_{g}(+\infty)$
$\displaystyle=\frac{\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}}{\Delta_{\text{rd}}-\mathrm{i}\kappa/2},$
(28) $\displaystyle\alpha_{e}(+\infty)$
$\displaystyle=\frac{\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}}{\Delta_{\text{rd}}+\chi_{\text{qr}}-\mathrm{i}\kappa/2},$
(29) $\displaystyle\alpha_{f}(+\infty)$
$\displaystyle=\frac{\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}}{\Delta_{\text{rd}}+2\chi_{\text{qr}}-\mathrm{i}\kappa/2},$
(30)

which can be anticipated from the solutions of the quantum Langevin equation
(QLE) when the resonator is driven by a classical source [30] , i.e.,

$\dot{\hat{a}}(t)=-\mathrm{i}\Big{[}\Delta_{\text{rd}}+\big{(}\chi_{\text{qr}}\ket{e}\\!\bra{e}+2\chi_{\text{qr}}\ket{f}\\!\bra{f}\big{)}\Big{]}\hat{a}(t)\\\
-\frac{\kappa}{2}\hat{a}(t)+\mathrm{i}\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}.$
(33)

What is not obvious by looking at the QLE is the dephasing rate captured in
Eq.(25)-(27). Note that the solution is obtained in the rotating frame; to go
back to the rest frame, we just need to restore the phase
$e^{-\mathrm{i}\omega_{\text{d}}t}$ in each coherent state.

We can go one step further by tracing out the resonator part; in other words,
the reduced density operator for the qutrit is given by

$\displaystyle\hat{\rho}_{\mathcal{S}}(t)$
$\displaystyle=\Tr_{\mathcal{R}}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\rho}_{\mathcal{SR}}(t)\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}$
$\displaystyle=\sum_{a}p_{a}(0)\ket{a}\\!\bra{a}+\sum_{a,b}c_{ab}(t)\ket{a}\\!\bra{b};$
(34)

hence, $c_{ab}$ are simply the coherence of the qutrit (i.e., the off-diagonal
terms of the reduced density operator of the qutrit) and can be solved from
Eq.(22)-(27).

Before discussing the effective qutrit master, we briefly mention the case
when the thermal bath is equilibrated at a nonzero temperature. The operator
differential equations of $\hat{\rho}_{ab}$ are almost the same as before
except that we replace
$\kappa\mathcalboondox{D}\big{[}\hat{a}\big{]}\hat{\rho}_{ab}$ with
$\kappa\big{(}\bar{N}+1\big{)}\mathcalboondox{D}\big{[}\hat{a}\big{]}\hat{\rho}_{ab}+\kappa\bar{N}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{a}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{ab}$
for $\bar{N}>0$. Consequently, each scalar differential equations of $P_{ab}$
acquires a second partial derivative, i.e.,

$\dot{P}_{ab}=\big{(}\text{terms from the case $\bar{N}=0$, see Appendix
\ref{A.1}}\big{)}\\\
+\kappa\bar{N}\frac{\partial^{2}}{\partial\alpha\partial\beta}P_{ab}.$ (37)

Since Eq.(37) with $a=b$ has the same form as the classical Fokker-Planck
equation [26, 31], its solution is described by a 2D Gaussian distribution of
finite width. In other words, instead of building up a coherent state (which
is a delta function in the positive $P$-representation) in the resonator, the
external drive will excite a Gaussian state with a quadrature uncertainty
broadened by the thermal bath. This also means that the resonator state is now
a linear combination of a continuum of coherent states with amplitudes near
$\alpha_{g}$, $\alpha_{e}$, or $\alpha_{f}$. In contrast, if the bath is in
the vacuum state, a coherent state excited in the resonator will remain
coherent indefinitely. See Appendix A.2 for more details.

###  III.3 Effective Qutrit/Qudit Master Equation

Going back to the zero-temperature assumption, the fact that Eq.(III.2)
describes the ensemble-averaged time evolution of a qutrit suggests that we
can construct an effective qutrit master equation for dispersive measurement.
To relax the long-$T_{1}$ assumption, we adopt the method of displacement
operator used in [12] to solve for the qutrit case.

We move to the displaced frame by introducing a qutrit-state-dependent
displacement operator

$\hat{\mathsf{P}}(t)=\hat{\Pi}_{g}\hat{D}(\alpha_{g}(t))+\hat{\Pi}_{e}\hat{D}(\alpha_{e}(t))+\hat{\Pi}_{f}\hat{D}(\alpha_{f}(t)),$
(38)

where $\hat{\Pi}_{a}=\ket{a}\\!\bra{a}$ is the projection operator onto the
qutrit state $\ket{a}$ for $a\in\\{g,e,f\\}$ and

$\hat{D}(\alpha_{a}(t))=\exp[\alpha_{a}(t)\hat{a}_{\text{r}}^{\dagger}-\alpha_{a}^{*}(t)\hat{a}_{\text{r}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}$
(39)

is the (time-dependent) resonator displacement operator which displaces a
coherent state $\ket{\alpha_{a}(t)}$ to the vacuum state [27]. In other words,
by performing the transformation
$\hat{\rho}_{\mathcal{SR}}^{\mathsf{P}}=\mathsf{P}^{\dagger}\hat{\rho}_{\mathcal{SR}}\mathsf{P}$
on the combined state and
$\hat{O}^{\mathsf{P}}=\mathsf{P}^{\dagger}\hat{O}^{\mathsf{P}}\mathsf{P}$ for
any operator $\hat{O}$ in the laboratory frame, we have effectively removed
all the resonator photons entangled with the qutrit.

Let $\rho_{\mathcal{S},ab}=\bra{a}\hat{\rho}_{\mathcal{S}}\ket{b}$ with
$a,b\in\\{g,e,f\\}$ be the matrix element of the qutrit density operator in
the laboratory frame. As shown in Appendix B and C, by first solving the
matrix elements of the density operator in the displaced frame and then moving
back to the laboratory frame with the resonator state traced out, one arrives
at

$\displaystyle\dot{\rho}_{\mathcal{S},gg}$
$\displaystyle=\gamma_{1,ge}\rho_{\mathcal{S},ee}+\gamma_{1,gf}\rho_{\mathcal{S},ff},$
(40) $\displaystyle\dot{\rho}_{\mathcal{S},ee}$
$\displaystyle=-\gamma_{1,ge}\rho_{\mathcal{S},ee}+\gamma_{1,ef}\rho_{\mathcal{S},ff},$
(41) $\displaystyle\dot{\rho}_{\mathcal{S},ff}$
$\displaystyle=-(\gamma_{1,gf}+\gamma_{1,ef})\rho_{\mathcal{S},ff}.$ (42)
$\displaystyle\dot{\rho}_{\mathcal{S},ge}$
$\displaystyle=\Big{[}\mathrm{i}\bar{\omega}_{eg}-\gamma_{1,ge}/2-\gamma_{\phi,ge}+\Gamma_{\text{d},ge}\Big{]}\rho_{\mathcal{S},ge},$
(43) $\displaystyle\dot{\rho}_{\mathcal{S},gf}$
$\displaystyle=\Big{[}\mathrm{i}\bar{\omega}_{gf}-\gamma_{1,ge}/2-\gamma_{\phi,gf}+\Gamma_{\text{d},gf}\Big{]}\rho_{\mathcal{S},gf},$
(44) $\displaystyle\dot{\rho}_{\mathcal{S},ef}$
$\displaystyle=\Big{[}\mathrm{i}\bar{\omega}_{ef}-\gamma_{1,ef}/2-\gamma_{\phi,ef}+\Gamma_{\text{d},ef}\Big{]}\rho_{\mathcal{S},ef},$
(45)

where the extra dephasing rates appearing in Eq.(43)-(45) are given by

$\displaystyle\Gamma_{\text{d},ge}(t)=\Gamma_{\text{d},eg}(t)$
$\displaystyle=\chi_{\text{qr}}\imaginary(\alpha_{g}\alpha_{e}^{*}),$ (46)
$\displaystyle\Gamma_{\text{d},gf}(t)=\Gamma_{\text{d},fg}(t)$
$\displaystyle=2\chi_{\text{qr}}\imaginary(\alpha_{g}\alpha_{f}^{*}),$ (47)
$\displaystyle\Gamma_{\text{d},ef}(t)=\Gamma_{\text{d},fe}(t)$
$\displaystyle=\chi_{\text{qr}}\imaginary(\alpha_{e}\alpha_{f}^{*}),$ (48)

and the time evolution of $\alpha_{a}(t)$ are still govened by Eq.(22)-(24).
In addition,
$\bar{\omega}_{ba}=\tilde{\omega}_{b}-\tilde{\omega}_{a}-\Delta_{\text{d},ba}$
are the new transition frequencies with the extra shifts

$\displaystyle\Delta_{\text{d},eg}(t)=-\Delta_{\text{d},ge}(t)$
$\displaystyle=\chi_{\text{qr}}\real(\alpha_{e}\alpha_{g}^{*}),$ (49)
$\displaystyle\Delta_{\text{d},fg}(t)=-\Delta_{\text{d},gf}(t)$
$\displaystyle=2\chi_{\text{qr}}\real(\alpha_{f}\alpha_{g}^{*}),$ (50)
$\displaystyle\Delta_{\text{d},fe}(t)=-\Delta_{\text{d},ef}(t)$
$\displaystyle=\chi_{\text{qr}}\real(\alpha_{f}\alpha_{e}^{*}).$ (51)

The dephasing rate $\Gamma_{\text{d},ab}$ and the frequency shift
$\Delta_{\text{d},ba}$ are nothing more than the real and imaginary parts of
the last terms in Eq.(25)-(27).

On the one hand, since we have included the effect of $T_{1}$, Eq.(40)-(42)
simply restate the semiclassical rate equations. On the other hand, the time
evolution of the coherence shown in Eq.(43)-(45) is exactly the same as that
in Eq.(25)-(27). $\Gamma_{\text{d},ab}$ is known as the measurement-induced
dephasing [10] and is a function of the dispersive shift $\chi_{\text{qr}}$,
resonator decay rate ${\kappa}$, readout detuning $\Delta_{\text{rd}}$, and
readout amplitude $\epsilon$. Furthermore, the same conclusion holds for a
general qudit measured dispersively: When subject to a coherent readout probe,
the resonator is excited to a superposition of $D$ coherent states each
entangled with a qudit eigenstate. The complex amplitude of the $j$th (with
$j=0,...,D-1$) coherent state evolves according to

$\dot{\alpha}_{j}=-\mathrm{i}(\Delta_{\text{rd}}+\chi_{j}-\mathrm{i}\kappa/2)\alpha_{j}+\mathrm{i}\epsilon;$
(52)

that is, the readout probe sees a resonator with a frequency shift $\chi_{j}$
relative to the bare frequency $\omega_{\text{r}}$. In addition, the time
evolution of the qutrit populations, again, follows the rate equation, i.e.,

$\dot{\rho}_{\mathcal{S},jj}=\sum_{k>j}\gamma_{1,jk}\rho_{\mathcal{S},kk}-\sum_{k<j}\gamma_{1,jk}\rho_{\mathcal{S},jj},$
(53)

where $\rho_{\mathcal{S},jj}$ for $j=0,...,D-1$ is the population of the $j$th
qudit eigenstate and $\gamma_{1,jk}$ with $j<k$ is the decay rate from state
$\ket{k}$ to $\ket{j}$. The time evolution of coherence term
$\rho_{\mathcal{S},jk}$ ($j<k$) is given by

$\dot{\rho}_{\mathcal{S},jk}=\Big{[}\mathrm{i}\bar{\omega}_{kj}-\gamma_{1,jk}/2-\gamma_{\phi,jk}\Gamma_{\text{d},jk}\Big{]}\rho_{\mathcal{S},jk},$
(54)

where $\gamma_{\phi,jk}$ is the pure dephasing,

$\Gamma_{\text{d},jk}(t)=\Gamma_{\text{d},kj}(t)=(\chi_{k}-\chi_{j})\imaginary(\alpha_{j}\alpha_{k}^{*})$
(55)

the measurement-induced dephasing, and

$\bar{\omega}_{kj}(t)=-\bar{\omega}_{jk}(t)=(\tilde{\omega}_{k}-\tilde{\omega}_{j})+(\chi_{k}-\chi_{j})\real(\alpha_{k}\alpha_{j}^{*})$
(56)

the shifted transition frequencies.

At this point, one might attempt to write down an effective master equation
for the qutrit based on Eq.(40)-(42) and (43)-(45); however, the qutrit
frequencies appeared in Eq.(43)-(45), in general, do not satisfy the relation

$\bar{\omega}_{fg}-\bar{\omega}_{fe}=\bar{\omega}_{ge}$ (57)

for a three-level system, so we cannot write down an exact master equation of
the Lindblad form. Such a problem does not appear in the qubit case [12] since
the single transition frequency $\bar{\omega}_{ge}$ is not subject to any
constraint. Nevertheless, if the frequency shifts are much smaller than other
rate parameters, we can still approximate the qutrit alone as a simple
Markovian system, thus writing down an effective master equation

$\displaystyle\dot{\hat{\rho}}_{\mathcal{S}}$
$\displaystyle=-\frac{\mathrm{i}}{\hbar}\Big{[}\hat{H}_{\text{q,eff}},\hat{\rho}_{\mathcal{S}}\Big{]}+\gamma_{1,ge}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ge}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}}+\gamma_{1,gf}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{gf}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}}+\gamma_{1,ef}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ef}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}}$
$\displaystyle\ \ \ \
+\frac{\gamma_{\phi,ge}+\Gamma_{\text{d},ge}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ge}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}}+\frac{\gamma_{\phi,gf}+\Gamma_{\text{d},gf}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,gf}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}}+\frac{\gamma_{\phi,ef}+\Gamma_{\text{d},ef}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ef}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}},$
(58)

where the effective qutrit Hamiltonian is assumed to describe a self-
consistent set of energy levels (e.g., by ignoring the measurement-induced
frequency shifts). Moreover, for the same reason, it should be clear that we
cannot write down an exact effective master equation for any qudit with $D\geq
3$. However, by ignoring the frequency shifts, the effective master equation
for a qudit will take the form

$\displaystyle\dot{\hat{\rho}}_{\mathcal{S}}$
$\displaystyle=-\frac{\mathrm{i}}{\hbar}\Big{[}\hat{H}_{\text{q,eff}},\hat{\rho}_{\mathcal{S}}\Big{]}+\sum_{j=0}^{D-1}\sum_{k>j}\gamma_{1,jk}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{jk}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}}$
$\displaystyle\ \ \ \ \ \ \
+\sum_{j=0}^{D-1}\sum_{k>j}\frac{\gamma_{\phi,jk}+\Gamma_{\text{d},jk}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,jk}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{S}}.$
(59)

###  III.4 Measurement-Induced Dephasing

We have already pointed out the connection between the diagonal terms of the
effective master equation and the semi-classical rate equation. The more
interesting phenomenon lies in the time evolution of the coherence terms. In
particular, we see that the product of the dispersive shift and the imaginary
part of $\alpha_{a}\alpha_{b}^{*}$ induces a dephasing for each energy level
of the qutrit. There are three factors that affect the measurement-induced
dephasing rates:

  1. (i)

A large readout drive leads to large coherent state amplitudes and thus a
stronger dephasing during the measurement. From the point of view of
information theory, since the field leaks out from the resonator will have a
larger amplitude as well, we are more likely to gain useful information from
the measurement since the signal-to-noise ratio increases as the power of the
readout signal goes up 111However, remember that the dispersive coupling is
valid only in the low-photon case, so the readout power is usually not a good
parameter to adjust.. Nevertheless, our measurement of the coherent states
leaked out of the resonator is subject to the quadrature uncertainty, thus,
the measurement result will be distributed as Gaussians centered at
$\alpha_{g}$, $\alpha_{e}$, or $\alpha_{f}$. The uncertainty in the
measurement will then lead to a random backaction on the qutrit conditioned on
the measurement results. It is this random backaction that leads to the
dephasing of the qutrit.

  2. (ii)

Related to the coherent state amplitudes is the readout frequency. As shown in
Eq.(28)-(30), the field amplitudes are Lorentzian functions of the detuning.
Hence, for the same drive strength $\bar{a}_{\text{in}}$, the amplitude built
up inside the resonator will be the highest when the detuning is zero.
However, we cannot drive all three dressed frequencies with zero detuning
simultaneously, which means that we have to play with the readout frequency so
that the separations among the three coherent states are maximized for state
classification. Figure 1(b) shows the steady-state amplitudes of the resonator
for some arbitrary readout frequency near the bare resonator frequency; in
general, the complex steady-state amplitudes lie on a circle that goes through
the origin of the phase plane. In addition, Figure 1(c) plots the distance
between two coherent state amplitudes as a function of the readout frequency.
As shown in Section IV, the distance between two coherent states determines
the measurement rate.

  3. (iii)

A larger dispersive shift $\chi_{\text{qr}}$ will also lead to a faster
decoherence time. This, again, can be argued from an information-theory point
of view. The dispersive shift determines how well we can separate the three
qutrit-dependent resonator frequencies; hence the larger the dispersive
shifts, the easier the state classification. However, as we have seen in the
analysis of the qudit-resonator coupling, $\chi_{\text{qr}}$ is proportional
to the square of the coupling coefficient and is approximately inversely
proportional to the square of the detuning. For the dispersive coupling to be
valid, one cannot make the qutrit-resonator coupling arbitrarily large, thus
limiting the amount of $\chi_{\text{qr}}$ realizable in practice. Furthermore,
there is another ratio we can design to improve the state classification – the
ratio between the dispersive shift $\chi_{\text{qr}}$ and cavity decay rate
$\kappa$. The effect of $\kappa$ is hidden in the expression of the steady-
state amplitudes. As shown in Figure 1(d), the distance between the coherent
states can be improved by making $\chi_{\text{qr}}>\kappa$.

##  IV Effective Qutrit Stochastic Master Equation

An unconditioned master equation can be interpreted as the stochastic
trajectories ensemble-averaged over all the possible measurement outcomes
[33]. The combined system of the qudit and the resonator can be measured
either actively by us or implicitly by the environment. In Section III, we
ignored the information coming out of the resonator, which is equivalent to
assuming that all measurements are implicitly made by the environment. To
describe an active dispersive measurement by us on a qudit, we thus need to
retrieve the information that has been so far averaged out. Since measurements
are probabilistic in quantum mechanics, we introduce an effective SME to model
the random measurement outcomes.

###  IV.1 Heterodyne Detection in the Transmission Mode

Now, we consider a qutrit coupled to the cavity dispersively. To read out the
resonator state, we perform a heterodyne detection where the readout signal
coming out of the resonator travels on the transmission line (i.e., the
coaxial cables connecting the inside of the dilution refrigerator to the room-
temperature electronics) and is mixed at the room temperature with a strong
local oscillator (LO) signal

$\displaystyle\hat{V}_{\text{LO}}(t)$
$\displaystyle=\frac{\hat{a}_{\text{LO}}(t)e^{-\mathrm{i}\phi_{\text{LO}}}+\hat{a}^{\dagger}_{\text{LO}}(t)e^{\mathrm{i}\phi_{\text{LO}}}}{2}$
$\displaystyle\approx
V_{\text{LO}}\cos(\omega_{\text{LO}}t-\phi_{\text{LO}}),$ (60)

whose frequency $\omega_{\text{LO}}$ is different from the readout signal
$\omega_{\text{d}}$ by the intermediate frequency (IF) $\omega_{\text{IF}}$,
e.g., $\omega_{\text{IF}}=\omega_{\text{LO}}-\omega_{\text{d}}$. In fact, in a
typical IQ demodulation stage, the amplified readout signal is first divided
equally in power and then mixed separately by two LO signals whose phases are
$90^{\circ}$-out-of-phase. Subsequently, the analog IF signals, passing
through an analog-to-digital converter (ADC), are processed digitally and
demodulated finally to DC (zero frequency) as a complex number (so that one
can plot the measurement as a point in the phase plane). Unlike a homodyne
detection where $\omega_{\text{LO}}=\omega_{\text{d}}$, the heterodyne scheme
allows us to measure both quadratures of the field at the same time (but still
constrained by the uncertainty principle). In addition, by first moving to an
IF frequency ($\omega_{\text{IF}}\sim 100\text{ MHz}$ in our experiment), the
signal experiences less $1/f$ noise.

We will, however, not attempt to model the analog or digital demodulation part
of the heterodyne detection using quantum mechanics. Instead, we will work
directly with the coherent state coming out of the resonator and assume that
we can process the signal in the way described above and retrieve information
about the coherent states subject to quantum-mechanical noise and imperfect
measurement efficiency. For a fully quantum-mechanical description of the
output chain, including filtering and amplification, see [22].

Unlike the qubit measurement where the readout frequency is usually set to be
$\omega_{\text{r}}+\chi_{\text{qr}}/2$ so that the detunings seen by the
readout probe are exactly $\pm\chi_{\text{qr}}/2$ and the information can be
encoded in only one quadrature, there isn’t any symmetry we can utilize to
describe a general qudit measured using the heterodyne scheme. To analyze the
information encoded in the complex amplitude of a coherent state, we first
define (in the rotating frame of $\omega_{\text{d}}$)

$\hat{I}_{\phi}=\frac{\hat{a}e^{-\mathrm{i}\phi}+\hat{a}^{\dagger}e^{\mathrm{i}\phi}}{2}\
\ \text{ and }\ \
\hat{Q}_{\phi}=\frac{\hat{a}e^{-\mathrm{i}\phi}-\hat{a}^{\dagger}e^{\mathrm{i}\phi}}{2\mathrm{i}}$
(61)

to be the two quadrature operators of the resonator field, where $\phi$ models
the net phase coming from the cable delay and any rotation applied during data
processing. Similarly, for the coherent states introduced in Section III.2, we
define

$\displaystyle\bar{I}_{g}(t)$
$\displaystyle=\real(\alpha_{g}e^{-\mathrm{i}\phi}),\ \ \ \
\bar{Q}_{g}(t)=\imaginary(\alpha_{g}e^{-\mathrm{i}\phi}),$ (62)
$\displaystyle\bar{I}_{e}(t)$
$\displaystyle=\real(\alpha_{e}e^{-\mathrm{i}\phi}),\ \ \ \
\bar{Q}_{e}(t)=\imaginary(\alpha_{e}e^{-\mathrm{i}\phi}),$ (63)
$\displaystyle\bar{I}_{f}(t)$
$\displaystyle=\real(\alpha_{f}e^{-\mathrm{i}\phi}),\ \ \ \
\bar{Q}_{f}(t)=\imaginary(\alpha_{f}e^{-\mathrm{i}\phi}).$ (64)

We, again, first restricted to the qutrit case since the qudit case can be
generalized easily.

Note, however, the quadrature fields defined above live in the resonator and
what we observe is only the leakage of the resonator to the transmission line.
We denote the input and output traveling-wave annihilation operators at port
$i=1,2$ as $\hat{a}_{\text{in,out},i}$. Recall that associated with the QLE of
the resonator at port 2 is the boundary condition [13, 34, 35]

$\displaystyle\hat{a}_{\text{out}}(t)$
$\displaystyle\doteq\hat{a}_{\text{out},2}(t)$
$\displaystyle=-\hat{a}_{\text{in},2}(t)+\sqrt{\kappa_{\text{out}}}\hat{a}(t)\approx\sqrt{\kappa_{\text{out}}}\hat{a}(t),$
(65)

where we have dropped $\hat{a}_{\text{in},2}$ by assuming that the incoming
signal at port 2 (i.e., the output port) is isolated by a well-designed
circulator/isolator stage and is not amplified at the output stage (i.e., the
HEMT is approximately unilateral). Since in the transmission mode, the
resonator is a two-port device, there should be another boundary condition for
port 1 (i.e., the input port, see Figure 1(a)); in fact, we have been
implicitly using it to define the drive
$\epsilon=\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}$ entering at port 1 of
the resonator. Information leaking out from port 1 will be factored into the
quantum efficiency $\eta$ to be discussed.

In the Schrödinger picture, the boundary condition at port 2 implies that the
transmitted signal is in a coherent state with the complex amplitude

$\alpha_{\text{out}}(t)=\sqrt{\kappa_{\text{out}}}\alpha_{a}(t)$ (66)

if the resonator is in the coherent state $\ket{\alpha_{a}(t)}$ for
$a\in\\{e,g,f\\}$. Moreover, since
$\hat{a}_{\text{out}}^{\dagger}\hat{a}_{\text{out}}$ is the outgoing photon
flux, the average number of photons leaving the resonator from port 2 within
an infinitesimally short time $\Delta t$ is given by

$\bar{n}(t)=\kappa_{\text{out}}\Delta t|\alpha_{a}(t)|^{2}.$ (67)

In terms of the outgoing quadrature fields on the transmission line, we have

$\displaystyle I_{\text{out}}(t)=\sqrt{\kappa_{\text{out}}\Delta
t}\bar{I}_{a}(t),$ (68) $\displaystyle
Q_{\text{out}}(t)=\sqrt{\kappa_{\text{out}}\Delta t}\bar{Q}_{a}(t),$ (69)

where $\bar{I}_{a}$ and $\bar{Q}_{a}$ are approximately constant over short
$\Delta t$. More precisely, it is the signal

$\hat{V}_{\text{out}}(t)=\frac{\hat{a}_{\text{out}}(t)e^{-(\mathrm{i}\omega_{\text{d}}t+\phi)}+\hat{a}_{\text{out}}^{\dagger}(t)e^{\mathrm{i}(\omega_{\text{d}}t+\phi)}}{2}$
(70)

that is mixed with $\hat{V}_{\text{LO}}(t)$ and its out-of-phase signal
$\hat{V}_{\text{LO}}(t-\pi/2\omega_{\text{LO}})$. Nevertheless, it’s not hard
to see that Eq.(68) and (69), up to a constant scaling due to cable loss and
amplification, are effectively what one would get after filtering the higher
sideband of the mixer output at $\omega_{\text{d}}+\omega_{\text{LO}}$ and
demodulating the filtered signal digitally to DC. For this reason, we will
treat $I_{\text{out}}$ and $Q_{\text{out}}$ as the measurement outcome
directly and omit calculation involving $\hat{V}_{\text{LO}}$. Any loss due to
the nonideality of the mixer will be included in the quantum efficiency.

In reality, however, measurements are not perfect and not all the information
encoded in the photon flux can be captured; thus, the effective photon number
we can measure is only

$\bar{n}_{\text{eff}}(t)=\eta\kappa\Delta t|\alpha_{a}(t)|^{2},$ (71)

where $\eta\in[0,1]$ is the measurement efficiency. Since
$\eta_{\text{r}}=\kappa_{\text{out}}/\kappa<1$, the efficiency $\eta$ is
naturally lowered by $\eta_{\text{r}}$. Note that $\eta_{\text{r}}$ also
contains the effect of $\kappa_{\text{int}}$ since photons lost internally are
inaccessible to the detector. In addition, by performing a heterodyne
detection, we automatically half the efficiency due to the power division in
an IQ mixer. Moreover, note that when $\eta=0$, we would gain zero information
about the system and can only talk about the behavior of the qutrit in an
averaged sense; thus, the SME to be constructed should reduce to the
unconditioned master equation and we will verify this point later.

Unlike the mean amplitude, the uncertainty/variance associated with a coherent
state traveling on the transmission line is fixed (each quadrature has a
variance of $1/4$), so the signal-to-noise ratio is proportional to $\Delta
t$; in other words, the photon shot noise can be effectively reduced by
increasing the measurement time. Concretely, suppose the resonator is in one
of the coherent states $\ket{\alpha_{a}}$ associated with an energy eigenstate
$\ket{a}$ of the qutrit. Then, the conditional probability density of
measuring a particular point $(I,Q)$ in the phase plane with an integration
time of $\Delta t$ given the qutrit state $\ket{a}$ is a two-dimensional
Gaussian

$\displaystyle f(I,Q|\hat{\rho}_{\mathcal{S}}(t)=\ket{a}\\!\bra{a})$
$\displaystyle\\!\\!\propto\exp[-\frac{1}{2}\frac{(I-\sqrt{\eta\kappa\Delta
t}\,\bar{I}_{a})^{2}+(Q-\sqrt{\eta\kappa\Delta t}\,\bar{Q}_{a})^{2}}{1/4}]$
(72)

with a variance $1/4$ in each quadrature. We have also assumed that the two
arms of the mixer output have the same conversion loss to use a single $\eta$;
in other words, the measurement is balanced.

More generally, if the qutrit is in a superposition state, then the cavity
state, after tracing out the qutrit state is given by

$\displaystyle\hat{\rho}_{\mathcal{R}}(t)$
$\displaystyle=\Tr_{\mathcal{S}}[\hat{\rho}_{\mathcal{SR}}(t)]$
$\displaystyle=\sum_{a\in\\{g,e,f\\}}\\!\\!\rho_{\mathcal{S},aa}(t)\ket{\alpha_{a}(t)}\\!\bra{\alpha_{a}(t)},$
(73)

as suggested by Eq.(21). Hence, the total conditional probability density of
measuring $(I,Q)$ in the phase plane is now given by

$f(I,Q|\hat{\rho}_{\mathcal{SR}}(t))\propto\sum_{a\in\\{g,e,f\\}}\rho_{\mathcal{S},aa}(t)\exp[-\frac{1}{2}\frac{(I-\sqrt{\eta\kappa\Delta
t}\,\bar{I}_{a})^{2}+(Q-\sqrt{\eta\kappa\Delta t}\,\bar{Q}_{a})^{2}}{1/4}].$
(74)

Consequently, the entanglement generated by the dispersive coupling will
project the qutrit state to a new state (possibly mixed) based on the
measurement outcome $(I,Q)$ after $\Delta t$.

By introducing the integration time $\Delta t$, we have discretized the
continuous measurement into time steps $t_{0},t_{1},t_{2},...$ with a step
size of $\Delta t$. We formalize measurement and the backaction induced by the
measurement output by introducing, at each time step $t_{k}$, a continuum of
POVM [36, 26]

$\Big{\\{}\hat{E}_{IQ}(t_{k})=\hat{K}_{IQ}^{\dagger}(t_{k})\hat{K}_{IQ}(t_{k})\
|\ (I,Q)\in\mathbf{R}^{2}\Big{\\}}$ (75)

for the qutrit (i.e., the resonator is traced out in this description) with
the Kraus operators

$\displaystyle\hat{K}_{IQ}(t_{k})$
$\displaystyle=\mathscr{N}_{k}\\!\\!\sum_{a\in\\{g,e,f\\}}\\!\\!\exp\\{\\!-\Big{[}I-\sqrt{\eta\kappa\Delta
t}\,\bar{I}_{a}(t_{k})\Big{]}^{2}$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ -\Big{[}Q-\sqrt{\eta\kappa\Delta
t}\,\bar{Q}_{a}(t_{k})\Big{]}^{2}\Bigg{\\}}\hat{\Pi}_{a}$ (76)

for any point $(I,Q)$ in the phase plane. The normalization constant
$\mathscr{N}_{k}$ can be found by imposing

$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}I\mathrm{d}Q\,\hat{E}_{IQ}(t_{k})=1,$
(77)

as required by the completeness of the POVM. Using the Kraus operators, the
probability density of measuring $(I,Q)$ is

$\Tr[\hat{\rho}_{\mathcal{S}}(t_{k})\hat{E}_{IQ}(t_{k})]=\Tr[\hat{K}_{IQ}(t_{k})\hat{\rho}_{\mathcal{S}}(t_{k})\hat{K}_{IQ}^{\dagger}(t_{k})],$
(78)

which, of course, must agree with Eq.(74). Furthermore, the post-measurement
state conditioned on the outcome $(I,Q)$ is

$\hat{\rho}_{\mathcal{S}}(t_{k+1})=\hat{\rho}_{\mathcal{S}}(t_{k}+\Delta
t)=\frac{\hat{K}_{IQ}(t_{k})\hat{\rho}_{\mathcal{S}}(t_{k})\hat{K}_{IQ}^{\dagger}(t_{k})}{\Tr[\hat{K}_{IQ}(t_{k})\hat{\rho}_{\mathcal{S}}(t_{k})\hat{K}_{IQ}^{\dagger}(t_{k})]}.$
(79)

We emphasize that $\hat{\rho}_{\mathcal{S}}(t_{k+1})$ is the conditional
reduced density operator and thus is not the same as the reduced density
operator used in the effective qutrit master equation before. Nevertheless,
once averaged over all the possible measurement history, we should reproduce
the unconditional density operator. Furthermore, it should be clear from
Eq.(79) that the series of quantum channels form a Markov chain, making the
entire mathematical formalism easier to deal with.

###  IV.2 Heuristic Derivation of the Qutrit Stochastic Master Equation

Based on Eq.(79), we look for a stochastic differential equation in the
diffusive limit [37], i.e., $\Delta t\rightarrow 0$. As shown in Figure 2, the
three Gaussian clusters in Eq.(74) merge together to form approximately a new
Gaussian distribution as $\Delta t$ becomes sufficiently small [38]. Appendix
D uses this approximation (which is true to the first order in $\Delta t$) to
reduce the Kraus operator to

$\displaystyle\hat{K}_{IQ}(t_{k})$
$\displaystyle\approx\tilde{\mathscr{N}}_{k}\exp{-\Big{[}I-\sqrt{\eta\kappa\Delta
t}\hat{L}_{I}(t_{k})\Big{]}^{2}}$ $\displaystyle\ \ \ \ \ \ \ \
\times\exp{-\Big{[}Q-\sqrt{\eta\kappa\Delta t}\hat{L}_{Q}(t_{k})\Big{]}^{2}},$
(80)

in the diffusive limit, where the new operators

$\displaystyle\hat{L}_{I}(t)=\bar{I}_{g}(t)\hat{\Pi}_{g}+\bar{I}_{e}(t)\hat{\Pi}_{e}+\bar{I}_{f}(t)\hat{\Pi}_{f},$
(81)
$\displaystyle\hat{L}_{Q}(t)=\bar{Q}_{g}(t)\hat{\Pi}_{g}+\bar{Q}_{e}(t)\hat{\Pi}_{e}+\bar{Q}_{f}(t)\hat{\Pi}_{f}.$
(82)

act as the Lindblad operators for measuring the $I$ and $Q$ quadrature,
respectively.

Figure 2: Illustration of a weak measurement in the diffusive limit. As the
integration time of each sample reduces, the three Gaussian distributions
merge together and can be approximated by a single Gaussian. The mean vector
of the approximated distribution is given by the centroid of the three
original mean vectors (i.e., $(\bar{I}_{a},\bar{Q}_{a})$) weighted by the
probability of the corresponding qutrit state $\ket{a}$.

By expanding Eq.(79) using Eq.(80), one can derive (see Appendix D) the
effective SME

$\displaystyle\mathrm{d}\hat{\rho}$
$\displaystyle=\bigg{(}-\frac{\mathrm{i}}{\hbar}\Big{[}\hat{H}_{\text{q,eff}},\hat{\rho}\Big{]}+\gamma_{1,ge}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ge}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\gamma_{1,gf}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{gf}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle\ \ \ \ \ \ \ \
+\gamma_{1,ef}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ef}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\frac{\gamma_{\phi,ge}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ge}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle\ \ \ \ \ \ \ \
+\frac{\gamma_{\phi,gf}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,gf}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\frac{\gamma_{\phi,ef}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ef}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\bigg{)}\mathrm{d}t$
$\displaystyle\ \ \ \
+\left(\kappa\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{I}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\kappa\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{Q}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\right)\mathrm{d}t$
$\displaystyle\ \ \ \
+\sqrt{\eta\kappa}\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{I}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}W_{I}+\sqrt{\eta\kappa}\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{Q}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}W_{Q}$
(83)

for the conditional reduced density operator of the qutrit, with the
heterodyne measurement outcomes (i.e., the complex signal demodulated to DC)
encoded in two classical stochastic differential equations

$\displaystyle
V_{I}(t)=\sqrt{\eta\kappa}\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}2\hat{L}_{I}(t)\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}+\xi_{I}(t),$
(84) $\displaystyle
V_{Q}(t)=\sqrt{\eta\kappa}\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}2\hat{L}_{Q}(t)\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}+\xi_{Q}(t),$
(85)

where $\langle\hat{c}\rangle=\Tr(\hat{\rho}\hat{c})$. The outcomes
$V_{I,Q}(t)$ are proportional to $I$ and $Q$ signals measured by the ADC in
the real experiment, but are rescaled to remove $\sqrt{\Delta t}$ in Eq.(68)
and (69). Note that the subscript $\mathcal{S}$ is dropped with the
understanding that $\hat{\rho}$ represents the conditional density operator of
the qutrit only.

In the SME, the measurement superoperator
$\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}$
associated with an operator $\hat{L}$ is defined via [39, 37]

$\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}=\hat{L}\hat{\rho}+\hat{\rho}\hat{L}^{\dagger}-\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}+\hat{L}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}\hat{\rho}.$
(86)

In addition, $W_{I,Q}$ are two independent classical Wiener processes and
$\xi_{I,Q}(t)=\dot{W}_{I,Q}(t)$ are classical white-noise signals satisfying

$\displaystyle\mathbb{E}[\xi_{I}(t)]=\mathbb{E}[\xi_{Q}(t)]=\mathbb{E}[\xi_{I}(t)\xi_{Q}(t^{\prime})]=0,$
(87)
$\displaystyle\mathbb{E}[\xi_{I}(t)\xi_{I}(t^{\prime})]=\mathbb{E}[\xi_{Q}(t)\xi_{Q}(t^{\prime})]=\delta(t-t^{\prime}).$
(88)

Formally, we should use Itô’s rule, i.e.,
$\mathrm{d}W_{I}^{2}=\mathrm{d}W_{Q}^{2}=\mathrm{d}t$ and
$\mathrm{d}W_{I}\mathrm{d}W_{Q}=0$ almost surely [40].

###  IV.3 Measurement Rates v.s. Measurement-Induced Dephasing Rates

At first glance, the measurement-induced dephasing rates defined in
Eq.(46)-(46) seem to be missed by the SME shown in Eq.(83). In fact,
$\kappa\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{I}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\kappa\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{Q}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
did not show up in the any of the unconditioned master equations introduced
before, at least not obvious in its current form. Nevertheless, one can show
that

$\displaystyle\kappa\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{I}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\kappa\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{Q}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle=\frac{\Gamma_{\text{m},ge}}{4}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ge}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\frac{\Gamma_{\text{m},gf}}{4}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,gf}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle\ \ \ \ \ \ \ \
+\frac{\Gamma_{\text{m},ef}}{4}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ef}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
(89)

by a simple rearrangement of the dissipation superoperators. In Eq.(89),

$\Gamma_{\text{m},ab}(t)=\kappa|\beta_{ab}(t)|^{2}$ (90)

is the measurement rate associated with the operator $\hat{\sigma}_{z,ab}$ and
$\beta_{ab}(t)=\alpha_{a}(t)-\alpha_{b}(t)$ represents the vector connecting
the two resonator coherent states in the phase plane. Intuitively, the greater
the separation between any two coherent states or the larger the resonator
decay rate, the easier the state classification becomes and, thus, the higher
the measurement rate.

Moreover, it turns out that the same rates $\Gamma_{\text{m},ab}$ also appear
as the dephasing rates in the unconditioned master equation of the combined
system in the displaced frame (see Appendix B); thus, Eq.(83) can be related
to the unconditioned master equation in a heuristic sense.

Figure 3: A comparison of the measurement dephasing rates derived in the
laboratory frame and the measurement rates found in the displaced frame. Since
we have assumed a Markovian system in the derivation of the SME, the
measurement rate and measurement-induced dephasing rates are not exactly the
same in the transient parts. However, when the system reaches its steady
state, the two rates are equivalent.

Even though we have reproduced something that appeared in the displaced frame,
we have still not derived the exact measurement-induced dephasing rates
$\Gamma_{\text{d},ab}$ that appear in the unconditioned master equation of the
qutrit. However, if we choose to use the steady-state values of
$\alpha_{g,e,f}$ listed in Eq.(28)-(30), we can show analytically that

$\Gamma_{\text{m},ab}(+\infty)=2\Gamma_{\text{d},ab}(+\infty)$ (91)

for $a\neq b$ and, thus, the heterodyne measurement indeed induces a dephasing
at rate $\Gamma_{\text{d},ab}$ between two energy levels of the qutrit at
steady state. The steady-state behavior was observed in the qubit case [12],
but our results imply the generality of such a relationship for a qutrit (and
a qudit). Figure 3 provides a detailed comparison between
$\Gamma_{\text{m},ge}$ and $2\Gamma_{\text{d},ge}$ before reaching the steady
state. Such a discrepancy can exist since we have been assuming a Markovian
system; from the derivation of the effective qutrit master equation, however,
we know that the qutrit alone is not really Markovian and the information loss
in the transient of the resonator evolution is not fully gained by the
heterodyne detection.

To capture the full measurement-induced dephasing using the quantum trajectory
approach, one must include the resonator as part of the SME with the
measurement operators acting on the resonator directly. Appendix D.2 provides
a derivation of the effective qutrit SME from the SME of the combined system,
showing that we can simply replace $\Gamma_{\text{m},ab}(t)$ in the
decoherence terms of Eq.(83) with $2\Gamma_{\text{d},ab}(t)$ as expected.
Consequently, we arrive at the corrected effective qutrit stochastic master
equation

$\displaystyle\mathrm{d}\hat{\rho}$
$\displaystyle=\bigg{(}-\frac{\mathrm{i}}{\hbar}\Big{[}\hat{H}_{\text{q,eff}},\hat{\rho}\Big{]}+\gamma_{1,ge}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ge}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle\ \ \ \ \ \ \ \
+\gamma_{1,gf}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{gf}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}+\gamma_{1,ef}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ef}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle\ \ \ \ \ \ \ \
+\frac{\gamma_{\phi,ge}+\Gamma_{\text{d},ge}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ge}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle\ \ \ \ \ \ \ \
+\frac{\gamma_{\phi,gf}+\Gamma_{\text{d},gf}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,gf}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}$
$\displaystyle\ \ \ \ \ \ \ \
+\frac{\gamma_{\phi,ef}+\Gamma_{\text{d},ef}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ef}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\bigg{)}\mathrm{d}t$
$\displaystyle\ \ \ \
+\sqrt{\eta\kappa}\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{I}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}W_{I}+\sqrt{\eta\kappa}\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{Q}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}W_{Q},$
(92)

where $\eta<1/2$ as a consequence of the heterodyne detection. Moreover, since
$\mathbb{E}(\mathrm{d}W_{I,Q})=0$, we indeed reproduce the effective qutrit
master equation described in Section III.3 (at steady state) by taking the
ensemble average of the SME.

###  IV.4 Generalization to the Measurement of a Qudit

The diffusive SME for the qutrit can be generalized in the qudit case. Besides
adding more decoherence channels, we need to redefine the Lindblad operators
for measurement:

$\hat{L}_{I}=\sum_{j=0}^{D-1}\bar{I}_{j}\hat{\Pi}_{j},\ \ \text{ and }\ \
\hat{L}_{Q}=\sum_{j=0}^{D-1}\bar{Q}_{j}\hat{\Pi}_{j},$ (93)

where the quadrature amplitudes of the $j$th coherent state are given by

$\bar{I}_{j}=\real(\alpha_{j}e^{-\mathrm{i}\phi})\ \ \text{ and }\ \
\bar{Q}_{j}=\imaginary(\alpha_{j}e^{-\mathrm{i}\phi})$ (94)

for $j=0,...,D-1$. Then, the effective qudit SME is given by

$\displaystyle\mathrm{d}\hat{\rho}$
$\displaystyle=-\frac{\mathrm{i}}{\hbar}\Big{[}\hat{H}_{\text{q,eff}},\hat{\rho}\Big{]}\mathrm{d}t+\sum_{j=0}^{D-1}\sum_{k>j}\gamma_{1,jk}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{jk}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}t$
$\displaystyle\ \ \ \ \ \ \
+\sum_{j=0}^{D-1}\sum_{k>j}\frac{\gamma_{\phi,jk}+\Gamma_{\text{d},jk}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,jk}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}t$
$\displaystyle\ \ \ \ \ \ \
+\sqrt{\eta\kappa}\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{I}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}W_{I}+\sqrt{\eta\kappa}\mathcalboondox{M}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{L}_{Q}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}\,\mathrm{d}W_{Q},$
(95)

where $\hat{\rho}$ represents the density operator of the qudit conditioned on
the measurement history. The measurement outcomes are still governed by
Eq.(84) and (85) but $\hat{L}_{I}$ and $\hat{L}_{Q}$ are now defined by
Eq.(93). Finally, taking the ensemble average of Eq.(95) gives the
unconditioned master equation stated in Eq.(59).

Figure 4: Monte Carlo simulation of the qutrit stochastic master equation with
an initial state given by Eq.(98). Three sample trajectories are shown in the
first three columns. Then, a thousand trajectories are averaged to produce the
last column. The first two rows plot the matrix elements of the qutrit density
operator. The last row shows the von Neumann entropy
$S(\hat{\rho})=-\Tr(\hat{\rho}\ln\hat{\rho})$, which peaks when the coherence
drops to zero. This is because the measurement-induced dephasing will first
erase the phase information in a superposition state; however, since the state
will be steered towards one of the energy eigenstates, the uncertainty drops
to zero eventually.

###  IV.5 Unraveling of the SME

In Eq.(92) and (95), the heterodyne measurement is modeled using stochastic
processes while the other decoherence terms are analyzed by taking the
ensemble average. Instead of having a mixture of quantum trajectories and
averaged time evolution in the same equation, one can unravel the SME to write
down a stochastic Schrödinger equation for the wave function. Since a mixed
state is not allowed in a Schrödinger equation, all the decoherence channels
must be described stochastically, i.e., by some random processes, resulting in
infinitely many realizations of the wave function. The unravelings of an SME
are not unique, but they must reproduce the decoherence terms in the SME after
averaging over all the realizations of the random processes associated with
the decoherence channels.

One possible unraveling uses the quantum state diffusion mode, which works
naturally with the diffusive limit used in describing the continuous
measurement. Given the set of pairs
$\mathopen{\scalerel*[5.5pt]{\Big{\\{}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\\{}}}}}\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{L}_{i},\gamma_{i}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}\mathclose{\scalerel*[5.5pt]{\Big{\\}}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\\}}}}}}_{i}$
with Lindblad operators $\hat{L}_{i}$ and the associated rates $\gamma_{i}$,
the quantum state diffusion equation is given by

$\displaystyle\mathrm{d}\ket{\Psi}$
$\displaystyle=-\frac{\mathrm{i}}{\hbar}\hat{H}_{\text{eff}}\ket{\Psi}\mathrm{d}t+\sum_{i}\sqrt{\gamma_{i}}\Big{(}\hat{L}_{i}-\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{i}\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}\Big{)}\mathrm{d}W_{i}$
$\displaystyle-\frac{1}{2}\sum_{i}\gamma_{i}\left(\hat{L}_{i}^{\dagger}\hat{L}_{i}+\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{i}\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{i}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}-2\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{i}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}\hat{L}_{i}\right)\ket{\Psi}\mathrm{d}t,$
(96)

where
$\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{i}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}=\bra{\Psi}\hat{L}_{i}^{\dagger}\ket{\Psi}$
is the expectation value computed for a single realization of $\ket{\Psi}$ and
$\\{W_{i}\\}_{i}$ are independent Wiener processes satisfying
$\mathbb{E}[\mathrm{d}W_{i}(t)]=\mathbb{E}[\mathrm{d}W_{i}(t)\mathrm{d}W_{j}(t)]=0$
for $i\neq j$ and $\mathbb{E}[\mathrm{d}W_{i}^{2}(t)]=\mathrm{d}t$. By
substituting

$\mathopen{\scalerel*[5.5pt]{\Big{\\{}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\\{}}}}}\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{L}_{i},\gamma_{i}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}\mathclose{\scalerel*[5.5pt]{\Big{\\}}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\\}}}}}}_{i}=\Big{\\{}\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{\sigma}_{jk},\gamma_{1,jk}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}},\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{\sigma}_{z,jk},\gamma_{\phi,jk}/2\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}\Big{\\}}_{k>j}$
(97)

into Eq.(96) and using the properties of the Wiener processes, one can easily
compute $\mathbb{E}[\mathrm{d}\hat{\rho}]$ and obtain the decoherence terms in
Eq.(95). In other words, the conditioning of the wave function on the
realizations of the Wiener processes is marginalized to create a mixed state
of the qudit.

Moreover, the random processes $W_{I}$ and $W_{Q}$ can be treated as a part of
Eq.(96) by adding
$\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{L}_{I},\sqrt{\eta\kappa}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}$,
$\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{L}_{I},\sqrt{(1-\eta)\kappa}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}$,
$\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{L}_{Q},\sqrt{\eta\kappa}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}$
and
$\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{L}_{Q},\sqrt{(1-\eta)\kappa}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}$
into Eq.(97) with independent Wiener processes $W_{I}$, $W_{I}^{\prime}$,
$W_{Q}$, and $W_{Q}^{\prime}$. Then, the full SME (in the steady state) can be
found by only taking the expectation with respect to $W_{I}^{\prime}$ and
$W_{Q}^{\prime}$, thus producing the correct measurement efficiency (see
Appendix D for modeling an imperfect measurement).

##  V Simulations and Experiments

###  V.1 Quantum Trajectories and Ensemble Averages

Eq.(92) can be simulated using finite difference just like a normal
differential equation. However, since the measurement is random, each finite-
difference simulation of Eq.(92) must be accompanied by a specific realization
of $W_{I}$ and $W_{Q}$. In other words, Eq.(92) can only be simulated in the
Monte-Carlo sense. In addition, since a finite-difference simulation
discretizes the time into steps of size $\Delta t$, one need to draw $\Delta
W_{I}$ and $\Delta W_{Q}$ from two independent Gaussian distributions both
with mean zero and variance $\Delta t$.

For example, consider an arbitrary qutrit state

$\hat{\rho}_{\mathcal{S}}(0)=\begin{pmatrix}0.5&0.3&0.36\\\ 0.3&0.2&0.24\\\
0.36&0.24&0.3\end{pmatrix}.$ (98)

If we repeat the simulation with this initial state a thousand times, we will
obtain a thousand distinct quantum trajectories of the qutrit. Three sample
trajectories are shown in Figure 4, along with the corresponding von Neumann
entropy. The last column of the figure gives the sample-averaged state and
entropy. Indeed, the sample-averaged population and coherence agree with the
ensemble-averaged ones shown in Figure 1(e)-(f), verifying that the
unconditioned master equation is the expectation of the SME.

What is not obvious from the unconditioned master equation is the convergence
of the qutrit state to one of the energy eigenstates as shown in the first row
of Figure 4. Although the qutrit is measured continuously and weakly, each
infinitesimal measurement can push the qutrit to a new state; consequently,
due to the nature of the dispersive coupling, the qutrit will slowly converge
to the pointer states $\ket{g}$, $\ket{e}$, and $\ket{f}$. For the qubit SME
in the long-$T_{1}$ limit, one can construct a Lyapunov function to show the
convergence of the qubit state to either $\ket{g}$ and $\ket{e}$ when the
readout probe is sent at $\omega_{\text{r}}+\chi_{\text{qr}}/2$ [37]. However,
for a qutrit or a qudit, a Lyapunov function is not known to the best of the
author’s knowledge. Despite the fact that each trajectory converges to an
energy eigenstate, the sample-averaged population is still flat as shown in
the last column of 4, which is a manifestation of the QND nature of dispersive
measurement. Since the dispersive measurements on average do not modify the
populations of each qutrit levels, we can conclude that the measurement
outcomes on average are a faithful reproduction of the underlying qutrit
state.

However, when $\gamma_{1,ab}$ are finite, the convergence of the qutrit is not
exactly true even if the quantum trajectories appear to have converged to an
eigenstate in a short amount of time. We can gain more insight by making the
simulation longer, which effectively increases the measurement time. As shown
in Figure 5, if we set the measurement time to be longer than
$1/\gamma_{1,ab}$, we would observe the shift of populations in the phase
plane. In other words, the sample-averaged populations, i.e., $\rho_{gg}$,
$\rho_{ee}$, and $\rho_{ff}$, are no longer a constant as appeared in Figure 4
if we extend the simulation time. Hence, in practice, the QND nature of the
measurement is limited by the decay rates and we cannot simply improve the
signal-to-noise ratio by increasing the measurement time arbitrarily.

Furthermore, by looking at the time evolution of the qutrit populations, we
observe jumps when the measurement time is sufficiently long. Similar behavior
of the quantum trajectories was discussed for the qubit case in [12]; for a
general qudit, as long as $\gamma_{1,jk}$ is nonzero, similar jumps will
happen so that the qudit state can eventually reach the ground state as
required by the unconditioned master equation.

Figure 5: Effect of qutrit decay in a long heterodyne measurement. A thousand
random quantum trajectories have been simulated with a measurement time
$T=40\text{ us}$ and with the qutrit parameters
$1/\gamma_{1,ge}=1/\gamma_{1,ef}=35\text{ us}$ and $1/\gamma_{1,gf}=1\text{
ms}$. a. For each quantum trajectory, the heterodyne measurement outcomes
$V_{I}(t)$ and $V_{Q}(t)$ are averaged over the measurement time $T$ and are
plotted as a point in the $I-Q$ plane. Besides the three Gaussian clusters
expected from a QND measurement, we also observe streams of points connecting
the three clusters, which represent the leakage of population from the excited
states to the ground states. b. When the measurement time is on the order of
the decay time, the random process exhibits a jumping characteristic on the
large time scale.

###  V.2 Measurement Shot Noise

Figure 6: Comparison between the simulation and experiment. The total cavity
decay rate is measured to be $\kappa/2\pi=2.7$ MHz and the dispersive shift is
$\chi_{\text{qr}}/2\pi=0.6$ MHz. The same parameters are used in the
simulation with the measurement efficiency set to $\eta=0.04$ (the experiment
is performed without a quantum-limited amplifier). The quadrature signals are
collected by AlazarTech digitizer with a sampling frequency of $1$ GHz, which
can be treated at $1/\Delta t$ in the derivation of the diffusive SME. a.
Changing the measurement time $T$. b. Sweeping the readout frequency
$f_{\text{d}}$ (i.e., the cavity drive frequency).

To validate the qutrit SME, real dispersive measurement is performed on a
transmon qutrit with $\tilde{\omega}_{\text{q}}/2\pi=4.48$ GHz and
$|\alpha_{\text{q}}|/2\pi\approx 280$ GHz coupled to a 3D aluminum cavity with
$\omega_{\text{r}}/2\pi=6.7835$ GHz when the qutrit is in the ground state.
The qutrit decay rates have been characterized to be larger than $30$ us for
$T_{1,ge}$, $T_{1,gf}$, and $T_{1,ef}$ and the maximum readout time used in
the experiment is much shorter than the decay time to ensure the validity of
the QND measurement. In addition, the shortest dephasing time is measured to
be $T_{2,gf}=1/\gamma_{2,gf}=3$ us, which does not play a significant role
since the measurement-induced dephasing happens at a much higher rate as can
be seen in the second row of Figure 4.

In all the subsequent experiments, the qutrit is excited to an equal
superposition state by two consecutive $\pi$-pulses to create three equally
weighted clusters in the phase plane. We include two types of comparison
between the experiment and theory. The first one is on the measurement time,
showing the reduction of the measurement uncertainty as more information leaks
out of the readout cavity. The second comparison tests how accurately the
theory can predict the steady-state amplitudes of the coherent states as a
function of the readout detuning.

In a real dispersive measurement, what we have access to are the transmitted
signals $V_{I}(t)$ and $V_{Q}(t)$. By calculating the time averages of
$V_{I}(t)$ and $V_{Q}(t)$ (with possibly weighting functions to account for
the transient behavior), denoted by $\bar{V}_{I}$ and $\bar{V}_{Q}$, we obtain
a complex number $\bar{V}_{I}+\mathrm{i}\bar{V}_{Q}$, which can be plotted on
a phase plane. Repeating the experiment a large number of times generates a
scattering plot, whose distribution reflects the qutrit state. One can then
compare the simulated $\bar{V}_{I}+\mathrm{i}\bar{V}_{Q}$ with the one
measured in experiments to verify the correctness of the SME. In Figure 6(a),
we vary the measurement time up to $3$ us to verify the effect of the shot
noise due to the quadrature of the coherent states and the broadening due to
the measurement inefficiency. As suggested by Eq.(84) and (85), by increasing
the measurement time and computing a more accurate time average of $V_{I,Q}$,
one will have access to the value of
$\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{I,Q}(t)\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}$
up to small variations in the transient part of $\hat{L}_{I,Q}(t)$. Since
$\hat{L}_{I,Q}(t)$ is a superposition of the projection operators weighted by
the quadrature amplitudes of the resonator coherent state, the qutrit state
(which determines the probability amplitude of each energy level) directly
affects the chance that one of the projection operator $\hat{\Pi}_{a}$ is
picked in the frequentist’s point of view. Hence,
$\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{I,Q}(t)\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}$
encodes the information about the qutrit state and a longer measurement time
effectively increases signal-to-noise ratio for measuring
$\mathopen{\scalerel*[5.5pt]{\Big{\langle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\langle}}}}}\hat{L}_{I,Q}(t)\mathclose{\scalerel*[5.5pt]{\Big{\rangle}}{\ensurestackMath{\addstackgap[1.5pt]{\big{\rangle}}}}}$
as illustrated in Figure 6(a).

###  V.3 Readout Frequency

With the readout time fixed at $3$ us, we sweep the readout frequency and
observe the change in the amplitude and phase of the (time-averaged)
transmitted signal. A comparison between the simulated and measured phase
planes is shown in Figure 6(b). As predicted by Figure 1(b), the center of
each Gaussian cluster lies on a circle that goes through the center.

When the readout frequency is far away from any of the state-dependent
resonator frequencies (i.e., $\omega_{\text{r}}+j\chi_{\text{qr}}$ for
$j=0,1,2$), the three Gaussian clusters overlap on one another, making the
state classification impossible. However, this also means that the qutrit
states are well protected (both in terms of the population and coherence)
since the rate at which information can leak out of the readout cavity at the
steady state is proportional to the separation between the clusters (see
Eq.(90) and (91)).

For a transmon-type qutrit where the dispersive shifts are $0$,
$\chi_{\text{qr}}$ and $2\chi_{\text{qr}}$, we can set the readout frequency
to be $\omega_{\text{r}}+\chi_{\text{qr}}$ so that the three Gaussian clusters
will be symmetrically placed in the phase plane as shown in the last column of
6(b). However, for a general qudit with unequally-spaced dispersive shifts,
there isn’t any symmetry one can utilize, which makes the SME simulation of
great usefulness.

Given that the simulation and experiment can be well matched, one can now find
an optimal set of design and experimental parameters that corresponds to a
large separation between the clusters and thus gives the best state
classification. In addition, other parameters such as measurement efficiency
can be determined for each experimental setup; subsequently, the simulation
used to design the qudits can adopt the same measurement efficiency to predict
the outcomes of the heterodyne detection.

##  VI Conclusion

In this work, we have introduced and solved the unconditioned master equation
describing the dispersive coupling between a qudit and a resonator.
Consequently, the concept of measurement-induced dephasing and frequency
shifts are extended to a qutrit and a general qudit. Two approaches – the
positive $P$-representation and the displaced frame – employed in the qubit
analysis have been adopted to solve the unconditioned master equation for the
qudit for the first time. Unlike the qubit case where the shift of the
transition frequency is unconstrained, we observe a non-Markovian nature of
the dispersive measurement for qudit with $D\geq 3$. In the presence of
negligible frequency shifts due to the readout probe, we arrive at an
effective (unconditioned) master equation for a qudit.

Given the analysis of the unconditioned qudit state, we then study the
conditioned qudit state, where the information leaving the combined system
(i.e., the transmitted readout signal) is retrieved and processed in both the
analog and digital domain. In particular, we focus on heterodyne detection,
where both quadratures of the coherent signals traveling along the
transmission line can be measured with the efficiency halved. In the diffusive
limit where each measurement has an infinitesimal integration time, we can
model the continuous monitoring of the qudit by a stochastic master equation
(SME). In addition, the measurement outcomes of the two quadratures are
described by two classical stochastic differential equations with the same
Wiener processes that appeared in the SME. In this way, the quantum dynamics
of the qudit are related to the classical information about the coherent
signal, allowing us to depict the stochastic nature of quantum measurement.

Finally, we compare the simulation of the SME with real experiments on a
transmon qutrit coupled to a readout cavity. The accuracy of the model is
tested with a sweep of the measurement time and of the readout frequency. Both
comparisons demonstrate the preciseness of the SME in predicting the formation
of clusters in the phase plane. Since our model does not rely on the type of
qudit, it can be applied to the measurement of other new superconducting qudit
topologies, such as a multimon. In practice, the analytical solution to the
unconditioned master SMEs can be critical guidance to the qudit design and
readout, enabling rapid prediction of the location of the clusters in the
phase plane and providing a potential theoretical foundation for understanding
the role of weak measurement in quantum feedback control.

###### Acknowledgements.

All simulations of the density operators were performed using QuTiP and NumPy.
The work at the University of California, Los Angeles (UCLA) is supported by
the National Science Foundation (NSF) Graduate Research Fellowship Program
under grant number DGE-2034835, the NSF Graduate Research Traineeship on
Quantum Science and Engineering at UCLA under grant number DGE-2125924, the
Office of Naval Research, and the Army Research Office. Any opinions,
findings, and conclusions or recommendations expressed in this material are
those of the author(s) and do not necessarily reflect the views of the NSF,
the Office of Naval Research, and the Army Research Office. The work at
Lawrence Livermore National Laboratory is performed under Contract No. DE-
AC52-07NA27344 by the U.S. Department of Energy.

##  Appendix A Solving the Master Equation of the Combined System (Qutrit +
Resonator)

###  A.1 Zero Temperature

To solve the master equation, we project the density operator of the composite
system onto the energy eigenbasis of the qutrit and thus introduce the
operators

$\hat{\rho}_{ab}(t)=\bra{a}\hat{\rho}_{\mathcal{SR}}(t)\ket{b}$ (99)

for $a,b\in\\{g,e,f\\}$; in other words, the reduced density operator can be
decomposed into

$\displaystyle\hat{\rho}_{\mathcal{SR}}(t)$
$\displaystyle=\hat{\rho}_{gg}\ket{g}\\!\bra{g}+\hat{\rho}_{ge}\ket{g}\\!\bra{e}+\hat{\rho}_{gf}\ket{g}\\!\bra{f}$
$\displaystyle\ \ \ \
+\hat{\rho}_{eg}\ket{e}\\!\bra{g}+\hat{\rho}_{ee}\ket{e}\\!\bra{e}+\hat{\rho}_{ef}\ket{e}\\!\bra{f}$
$\displaystyle\ \ \ \
+\hat{\rho}_{fg}\ket{f}\\!\bra{g}+\hat{\rho}_{fe}\ket{f}\\!\bra{e}+\hat{\rho}_{ff}\ket{f}\\!\bra{f},$
(100)

and, after expanding the master equation using the nine operators, we obtain
nine coupled operator differential equations

$\displaystyle\dot{\hat{\rho}}_{gg}$
$\displaystyle=-\mathrm{i}\Delta_{\text{rd}}\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{gg}\right]$
$\displaystyle\ \ \ \
+\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{gg}\right]+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{gg},$
(101) $\displaystyle\dot{\hat{\rho}}_{ee}$
$\displaystyle=-\mathrm{i}(\chi_{\text{qr}}+\Delta_{\text{rd}})\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{ee}\right]$
$\displaystyle\ \ \ \
+\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{ee}\right]+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{ee},$
(102) $\displaystyle\dot{\hat{\rho}}_{ff}$
$\displaystyle=-\mathrm{i}(2\chi_{\text{qr}}+\Delta_{\text{rd}})\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{ff}\right]$
$\displaystyle\ \ \ \
+\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{ff}\right]+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{ff},$
(103) $\displaystyle\dot{\hat{\rho}}_{ge}$
$\displaystyle=\mathrm{i}\tilde{\omega}_{\text{q}}\hat{\rho}_{ge}-\mathrm{i}\Delta_{\text{rd}}\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{ge}\right]\\!+\mathrm{i}\chi_{\text{qr}}\hat{\rho}_{ge}\hat{a}^{\dagger}\hat{a}$
$\displaystyle\ \ \ \
+\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{ge}\right]\\!+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{ge}-\gamma_{2,ge}\hat{\rho}_{ge},$
(104) $\displaystyle\dot{\hat{\rho}}_{eg}$
$\displaystyle=-\mathrm{i}\tilde{\omega}_{\text{q}}\hat{\rho}_{eg}+\mathrm{i}\Delta_{\text{rd}}\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{eg}\right]\\!-\mathrm{i}\chi_{\text{qr}}\hat{\rho}_{eg}\hat{a}^{\dagger}\hat{a}$
$\displaystyle\ \ \ \
-\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{eg}\right]\\!+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{eg}-\gamma_{2,eg}\hat{\rho}_{eg},$
(105) $\displaystyle\dot{\hat{\rho}}_{gf}$
$\displaystyle=\mathrm{i}(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\hat{\rho}_{gf}-\mathrm{i}\Delta_{\text{rd}}\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{gf}\right]\\!$
$\displaystyle\ \ \ \
+\mathrm{i}2\chi_{\text{qr}}\hat{\rho}_{gf}\hat{a}^{\dagger}\hat{a}+\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{gf}\right]\\!$
$\displaystyle\ \ \ \
+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{gf}-\gamma_{2,gf}\hat{\rho}_{gf},$
(106) $\displaystyle\dot{\hat{\rho}}_{fg}$
$\displaystyle=-\mathrm{i}(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\hat{\rho}_{fg}+\mathrm{i}\Delta_{\text{rd}}\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{fg}\right]\\!$
$\displaystyle\ \ \ \
-\mathrm{i}2\chi_{\text{qr}}\hat{\rho}_{fg}\hat{a}^{\dagger}\hat{a}-\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{fg}\right]\\!$
$\displaystyle\ \ \ \
+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{fg}-\gamma_{2,fg}\hat{\rho}_{fg},$
(107) $\displaystyle\dot{\hat{\rho}}_{ef}$
$\displaystyle=\mathrm{i}(\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\hat{\rho}_{ef}-\mathrm{i}\Delta_{\text{rd}}\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{ef}\right]\\!$
$\displaystyle\ \ \ \
+\mathrm{i}\chi_{\text{qr}}(2\hat{\rho}_{ef}\hat{a}^{\dagger}\hat{a}-\hat{a}^{\dagger}\hat{a}\hat{\rho}_{ef})$
$\displaystyle\ \ \ \
+\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{ef}\right]\\!+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{ef}-\gamma_{2,ef}\hat{\rho}_{ef},$
(108) $\displaystyle\dot{\hat{\rho}}_{fe}$
$\displaystyle=-\mathrm{i}(\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\hat{\rho}_{fe}+\mathrm{i}\Delta_{\text{rd}}\left[\hat{a}^{\dagger}\hat{a},\hat{\rho}_{fe}\right]\\!$
$\displaystyle\ \ \ \
-\mathrm{i}\chi_{\text{qr}}(2\hat{\rho}_{fe}\hat{a}^{\dagger}\hat{a}-\hat{a}^{\dagger}\hat{a}\hat{\rho}_{fe})$
$\displaystyle\ \ \ \
-\mathrm{i}\left[\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a},\hat{\rho}_{fe}\right]\\!+\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{fe}-\gamma_{2,fe}\hat{\rho}_{fe}.$
(109)

Note that each operator $\hat{\rho}_{ab}$ lives in an infinite dimensional
space since $\mathscr{H}_{\mathcal{R}}$ is a Fock space. Nevertheless, it’s
possible to find a closed-form solution by invoking the positive
$P$-representation

$\hat{\rho}_{ab}(t)=\int\mathrm{d}^{2}\alpha\int\mathrm{d}^{2}\beta\frac{\ket{\alpha}\\!\bra{\beta^{*}}}{\innerproduct{\beta^{*}}{\alpha}}P_{ab}(\alpha,\beta,t).$
(110)

The reader can verify that the action of the creation and annihilation
operators in the operator space can be translated to some simple operations in
the positive $P$-representation:

$\displaystyle\hat{a}\hat{\rho}(t)$ $\displaystyle\ \ \longrightarrow\ \
\alpha P(\alpha,\beta,t),$ (111) $\displaystyle\hat{a}^{\dagger}\hat{\rho}(t)$
$\displaystyle\ \ \longrightarrow\ \
\left(\beta-\frac{\partial}{\partial\alpha}\right)P(\alpha,\beta,t),$ (112)
$\displaystyle\hat{\rho}(t)\hat{a}^{\dagger}$ $\displaystyle\ \
\longrightarrow\ \ \beta P(\alpha,\beta,t),$ (113)
$\displaystyle\hat{\rho}(t)\hat{a}$ $\displaystyle\ \ \longrightarrow\ \
\left(\alpha-\frac{\partial}{\partial\beta}\right)P(\alpha,\beta,t).$ (114)

As an example, let us use Eq.(111)-(114) to transform Eq.(A.1) into a scalar
equation:

$\displaystyle\dot{P}_{gg}$
$\displaystyle=-\mathrm{i}\Delta_{\text{rd}}\left[\left(\beta-\frac{\partial}{\partial\alpha}\right)(\alpha
P_{gg})-\left(\alpha-\frac{\partial}{\partial\beta}\right)(\beta
P_{gg})\right]$ $\displaystyle\ \ \ \
+\mathrm{i}\bigg{[}\epsilon\left(\beta-\frac{\partial}{\partial\alpha}\right)P_{gg}+\epsilon^{*}\alpha
P_{gg}$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ -\epsilon\beta
P_{gg}-\epsilon^{*}\left(\alpha-\frac{\partial}{\partial\beta}\right)P_{gg}\bigg{]}$
$\displaystyle\ \ \ \ +\kappa\bigg{[}\alpha\beta
P_{gg}-\frac{1}{2}\left(\beta-\frac{\partial}{\partial\alpha}\right)(\alpha
P_{gg})$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
-\frac{1}{2}\left(\alpha-\frac{\partial}{\partial\beta}\right)(\beta
P_{gg})\bigg{]}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[\left(-\mathrm{i}\epsilon+\mathrm{i}\Delta_{\text{rd}}\alpha+\frac{\kappa\alpha}{2}\right)P_{gg}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[\left(\mathrm{i}\epsilon^{*}-\mathrm{i}\Delta_{\text{rd}}\beta+\frac{\kappa\beta}{2}\right)P_{gg}\right],$
(115)

where $\partial_{t}P_{gg}(\alpha,\beta,t)$ is shorthanded as $\dot{P}_{gg}$.
By applying the same procedure to the other eight operators, we get, in total,
nine coupled scalar differential equations

$\displaystyle\\!\\!\\!\\!\dot{P}_{gg}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{gg}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{gg}\right],$
(116) $\displaystyle\\!\\!\\!\\!\dot{P}_{ee}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}\chi_{\text{qr}}\alpha+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{ee}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}\chi_{\text{qr}}\beta-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{ee}\right],$
(117) $\displaystyle\\!\\!\\!\\!\dot{P}_{ff}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}2\chi_{\text{qr}}\alpha+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{ff}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}2\chi_{\text{qr}}\beta-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{ff}\right],$
(118) $\displaystyle\\!\\!\\!\\!\dot{P}_{ge}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{ge}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}\chi_{\text{qr}}\beta-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{ge}\right]$
$\displaystyle\ \ \ \ +\mathrm{i}\chi_{\text{qr}}\alpha\beta
P_{ge}+\mathrm{i}\tilde{\omega}_{\text{q}}P_{ge}-\gamma_{2,ge}P_{ge},$ (119)
$\displaystyle\\!\\!\\!\\!\dot{P}_{eg}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}\chi_{\text{qr}}\alpha+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{eg}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{eg}\right]$
$\displaystyle\ \ \ \ -\mathrm{i}\chi_{\text{qr}}\alpha\beta
P_{eg}-\mathrm{i}\tilde{\omega}_{\text{q}}P_{eg}-\gamma_{2,ge}P_{eg},$ (120)
$\displaystyle\\!\\!\\!\\!\dot{P}_{gf}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{gf}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}2\chi_{\text{qr}}\beta-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{gf}\right]$
$\displaystyle\ \ \ \ +\mathrm{i}2\chi_{\text{qr}}\alpha\beta
P_{gf}+\mathrm{i}(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})P_{gf}-\gamma_{2,gf}P_{gf},$
(121) $\displaystyle\dot{P}_{fg}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}2\chi_{\text{qr}}\alpha+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{fg}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{fg}\right]$
$\displaystyle\ \ \ \ -\mathrm{i}2\chi_{\text{qr}}\alpha\beta
P_{fg}-\mathrm{i}(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})P_{fg}-\gamma_{2,gf}P_{fg},$
(122) $\displaystyle\dot{P}_{ef}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}\chi_{\text{qr}}\alpha+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{ef}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}2\chi_{\text{qr}}\beta-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{ef}\right]$
$\displaystyle\ \ \ \ +\mathrm{i}\chi_{\text{qr}}\alpha\beta
P_{ef}+\mathrm{i}(\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})P_{ef}-\gamma_{2,ef}P_{ef},$
(123) $\displaystyle\dot{P}_{fe}$
$\displaystyle=\frac{\partial}{\partial\alpha}\left[(-\mathrm{i}\epsilon+\mathrm{i}2\chi_{\text{qr}}\alpha+\mathrm{i}\Delta_{\text{rd}}\alpha+\kappa\alpha/2)P_{fe}\right]$
$\displaystyle\ \ \ \
+\frac{\partial}{\partial\beta}\left[(\mathrm{i}\epsilon^{*}-\mathrm{i}\chi_{\text{qr}}\beta-\mathrm{i}\Delta_{\text{rd}}\beta+\kappa\beta/2)P_{fe}\right]$
$\displaystyle\ \ \ \ -\mathrm{i}\chi_{\text{qr}}\alpha\beta
P_{fe}-\mathrm{i}(\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})P_{fe}-\gamma_{2,ef}P_{fe}.$
(124)

It should be noted that the differential equations of $P_{ab}$ are usually of
the type of a Fokker-Planck equation, which also includes the diffusive terms
(i.e., the second partial derivatives with respect to $\alpha$ and $\beta$).
However, since we have assumed that $\bar{N}(\omega_{\text{r}})=0$, there is
no terms of the form

$\hat{a}^{\dagger}\hat{\rho}_{ab}\hat{a}\ \longrightarrow\ \
\left(\beta-\frac{\partial}{\partial\alpha}\right)\left(\alpha-\frac{\partial}{\partial\beta}\right)P(\alpha,\beta,t).$
(125)

Even in the case where $\bar{N}>0$, the method of the positive
$P$-representation will still work but a sharp coherent state (see below)
inside the cavity will broaden itself diffusively in the phase plane.

Although looking complicated, the nine coupled equations admit simple
trajectories in the complex planes of $\alpha$ and $\beta$. We use the ansatze

$\displaystyle P_{gg}(\alpha,\beta,t)$
$\displaystyle=\delta^{(2)}(\alpha-\alpha_{g}(t))\delta^{(2)}(\beta-\alpha_{g}^{*}(t)),$
(126) $\displaystyle P_{ee}(\alpha,\beta,t)$
$\displaystyle=\delta^{(2)}(\alpha-\alpha_{e}(t))\delta^{(2)}(\beta-\alpha_{e}^{*}(t)),$
(127) $\displaystyle P_{ff}(\alpha,\beta,t)$
$\displaystyle=\delta^{(2)}(\alpha-\alpha_{f}(t))\delta^{(2)}(\beta-\alpha_{f}^{*}(t))$
(128)

for the diagonal terms and

$\displaystyle\\!\\!\\!\\!P_{ge}(\alpha,\beta,t)$
$\displaystyle=c_{ge}(t)\delta^{(2)}(\alpha-\alpha_{g}(t))\delta^{(2)}(\beta-\alpha_{e}^{*}(t)),$
(129) $\displaystyle\\!\\!\\!\\!P_{eg}(\alpha,\beta,t)$
$\displaystyle=c_{eg}(t)\delta^{(2)}(\alpha-\alpha_{e}(t))\delta^{(2)}(\beta-\alpha_{g}^{*}(t)),$
(130) $\displaystyle\\!\\!\\!\\!P_{gf}(\alpha,\beta,t)$
$\displaystyle=c_{gf}(t)\delta^{(2)}(\alpha-\alpha_{g}(t))\delta^{(2)}(\beta-\alpha_{f}^{*}(t)),$
(131) $\displaystyle\\!\\!\\!\\!P_{fg}(\alpha,\beta,t)$
$\displaystyle=c_{fg}(t)\delta^{(2)}(\alpha-\alpha_{f}(t))\delta^{(2)}(\beta-\alpha_{g}^{*}(t)),$
(132) $\displaystyle\\!\\!\\!\\!P_{ef}(\alpha,\beta,t)$
$\displaystyle=c_{ef}(t)\delta^{(2)}(\alpha-\alpha_{e}(t))\delta^{(2)}(\beta-\alpha_{f}^{*}(t)),$
(133) $\displaystyle\\!\\!\\!\\!P_{fe}(\alpha,\beta,t)$
$\displaystyle=c_{fe}(t)\delta^{(2)}(\alpha-\alpha_{f}(t))\delta^{(2)}(\beta-\alpha_{e}^{*}(t))$
(134)

for off-diagonal terms. Each diagonal term $P_{aa}$ represents a single
coherent state whose amplitudes are specified by the two delta functions.
Plugging the ansatze into Eq.(A.1)-(A.1), we obtain the time evolution of the
coherent states

$\displaystyle\dot{\alpha}_{g}$
$\displaystyle=-\mathrm{i}(\Delta_{\text{rd}}-\mathrm{i}\kappa/2)\alpha_{g}+\mathrm{i}\epsilon,$
(135) $\displaystyle\dot{\alpha}_{e}$
$\displaystyle=-\mathrm{i}(\Delta_{\text{rd}}+\chi_{\text{qr}}-\mathrm{i}\kappa/2)\alpha_{e}+\mathrm{i}\epsilon,$
(136) $\displaystyle\dot{\alpha}_{f}$
$\displaystyle=-\mathrm{i}(\Delta_{\text{rd}}+2\chi_{\text{qr}}-\mathrm{i}\kappa/2)\alpha_{f}+\mathrm{i}\epsilon.$
(137)

Besides the time evolution brought by the coherent states, each non-diagonal
term $P_{ab}$ ($a\neq b$) is modulated by an envelope function $c_{ab}$. By
substituting the ansatze together with Eq.(135)-(137) into Eq.(119)-(124), we
deduce that

$\displaystyle\dot{c}_{ge}$
$\displaystyle=\mathrm{i}(\tilde{\omega}_{\text{q}}+\mathrm{i}\gamma_{2,ge})c_{ge}+\mathrm{i}\chi_{\text{qr}}\alpha_{g}\alpha_{e}^{*}c_{ge},$
(138) $\displaystyle c_{eg}$ $\displaystyle=c_{ge}^{*},$ (139)
$\displaystyle\dot{c}_{gf}$
$\displaystyle=\mathrm{i}(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}}+\mathrm{i}\gamma_{2,gf})c_{gf}+\mathrm{i}2\chi_{\text{qr}}\alpha_{g}\alpha_{f}^{*}c_{gf},$
(140) $\displaystyle c_{fg}$ $\displaystyle=c_{gf}^{*},$ (141)
$\displaystyle\dot{c}_{ef}$
$\displaystyle=\mathrm{i}(\tilde{\omega}_{\text{q}}+\alpha_{\text{q}}+\mathrm{i}\gamma_{2,ef})c_{ef}+\mathrm{i}\chi_{\text{qr}}\alpha_{e}\alpha_{f}^{*}c_{ef},$
(142) $\displaystyle c_{fe}$ $\displaystyle=c_{ef}^{*}.$ (143)

Given the arbitrary initial conditions, the detailed time evolution of
$\alpha_{a}$ and $c_{ab}$ can be solved numerically. After integrating the
delta functions in the complex plane of $\alpha$ and $\beta$, we arrive at the
general solution (in the rotating frame)

$\displaystyle\hat{\rho}_{\mathcal{SR}}(t)=\sum_{a\in\\{g,e,f\\}}p_{a}(0)\ket{a}\\!\bra{a}\otimes\ket{\alpha_{a}(t)}\\!\bra{\alpha_{a}(t)}$
$\displaystyle+\sum_{a,b\in\\{g,e,f\\}}\frac{c_{ab}(t)}{\bra{\alpha_{b}(t)}\ket{\alpha_{a}(t)}}\ket{a}\\!\bra{b}\otimes\ket{\alpha_{a}(t)}\\!\bra{\alpha_{b}(t)},$
(144)

where $p_{g,e,f}(0)$ are the initial populations in each energy eigenstate of
the qutrit. The steady-state (complex) amplitudes of the cavity coherent
states are given by

$\displaystyle\alpha_{g}(+\infty)$
$\displaystyle=\frac{\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}}{\Delta_{\text{rd}}-\mathrm{i}\kappa/2},$
(145) $\displaystyle\alpha_{e}(+\infty)$
$\displaystyle=\frac{\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}}{\Delta_{\text{rd}}+\chi_{\text{qr}}-\mathrm{i}\kappa/2},$
(146) $\displaystyle\alpha_{f}(+\infty)$
$\displaystyle=\frac{\sqrt{\kappa_{\text{in}}}\bar{a}_{\text{in}}}{\Delta_{\text{rd}}+2\chi_{\text{qr}}-\mathrm{i}\kappa/2}.$
(147)

In the Appendix B and C, we will approach the same problem with a different
technique; nevertheless, the result will be identical, except that we will
include $\gamma_{1,ab}$ for completeness.

###  A.2 Nonzero Temperature

When $\bar{N}>0$, the master equation takes the form

$\displaystyle\dot{\hat{\rho}}_{\mathcal{SR}}$
$\displaystyle=-\frac{\mathrm{i}}{\hbar}\left[\hat{H}_{\text{eff}},\hat{\rho}_{\mathcal{SR}}\right]+\kappa\big{(}\bar{N}+1\big{)}\mathcalboondox{D}[\hat{a}]\hat{\rho}_{\mathcal{SR}}$
$\displaystyle\ \ \ \
+\kappa\bar{N}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{a}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{\mathcal{SR}}+\frac{\gamma_{2,ge}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,ge}\big{]}\hat{\rho}_{\mathcal{SR}}$
$\displaystyle\ \ \ \
+\frac{\gamma_{2,gf}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,gf}\big{]}\hat{\rho}_{\mathcal{SR}}+\frac{\gamma_{2,ef}}{2}\mathcalboondox{D}\big{[}\hat{\sigma}_{z,ef}\big{]}\hat{\rho}_{\mathcal{SR}},$
(148)

in the long-$T_{1}$ limit. The operator differential equations of
$\hat{\rho}_{ab}$ are almost the same as before except that we replace
$\kappa\mathcalboondox{D}[\hat{a}]\hat{\rho}_{ab}$ with

$\kappa\big{(}\bar{N}+1\big{)}\mathcalboondox{D}[\hat{a}]\hat{\rho}_{ab}+\kappa\bar{N}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{a}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}_{ab}.$
(149)

Consequently, the scalar differential equations for the positive
$P$-representations $P_{ab}$ acquire the second partial derivatives mention
before, i.e.,

$\dot{P}_{ab}=\big{(}\text{terms from the case {\hbox{\bar{N} =
0}}}\big{)}+\kappa\bar{N}\frac{\partial^{2}}{\partial\alpha\partial\beta}P_{ab}.$
(150)

Since Eq.(150) with $a=b$ has the same form as the classical Fokker-Planck
equation, we use Gaussian distributions now as the new ansatze

$\displaystyle P_{gg}(\alpha,\beta,t)$ $\displaystyle=\frac{1}{\pi
N(t)}\exp\left\\{-\frac{1}{N(t)}\big{[}\alpha-\alpha_{g}(t)\big{]}\big{[}\beta-\alpha_{g}^{*}(t)\big{]}\right\\},$
(151) $\displaystyle P_{ee}(\alpha,\beta,t)$ $\displaystyle=\frac{1}{\pi
N(t)}\exp\left\\{-\frac{1}{N(t)}\big{[}\alpha-\alpha_{e}(t)\big{]}\big{[}\beta-\alpha_{e}^{*}(t)\big{]}\right\\},$
(152) $\displaystyle P_{ff}(\alpha,\beta,t)$ $\displaystyle=\frac{1}{\pi
N(t)}\exp\left\\{-\frac{1}{N(t)}\big{[}\alpha-\alpha_{f}(t)\big{]}\big{[}\beta-\alpha_{f}^{*}(t)\big{]}\right\\},$
(153)

where we require that $\alpha_{g,e,f}$ still satisfy Eq.(135), (136), and
(137), respectively. Substituting the Gaussian distributions into the Fokker-
Planck equation, we obtain the differential equation of the variance $N$ (more
precisely, the variance of the two-dimensional Gaussian is $N/2$)

$\dot{N}(t)=-\kappa\big{[}N(t)-\bar{N}\big{]}.$ (154)

Suppose the composite system was in thermal equilibrium with the bath before
receiving the drive $\varepsilon_{\text{d}}(t)$, then we will simply use
$N(+\infty)=\bar{N}$ in $P_{gg}$, $P_{ee}$, and $P_{ff}$. In other words,
instead of building up a coherent state in the resonator, the external drive
will excite a Gaussian state with a quadrature uncertainty broadened by the
thermal bath. This also means that the resonator is a continuum combination of
coherent states with amplitudes near $\alpha_{g}$, $\alpha_{e}$, and
$\alpha_{f}$. In contrast, if the bath is in a vacuum state, $N(+\infty)=0$
and a coherent state initially excited in the resonator will remain coherent.

Unlike the case where $\bar{N}=0$, the coherence $c_{ab}$ also depends on
$\alpha$ and $\beta$ now. This is because $c_{ab}$ are affected by an infinite
collection of coherent states around $\alpha_{g}$, $\alpha_{e}$, and
$\alpha_{f}$. Nevertheless, we expect $c_{ab}$ to vanish on similar timescales
set by Eq.(138)-(143).

##  Appendix B Deriving the Master Equation of the Combined System (Qutrit +
Resonator) in the Displaced Frame

###  B.1 Qutrit-State-Dependent Displacement

Due to dispersive coupling, each eigenstate of the qutrit is entangled with a
coherent state of the resonator. In the last section, we have found three
differential equations for the complex amplitudes $\alpha_{g}$, $\alpha_{e}$,
and $\alpha_{f}$ of the coherent states. To make this entanglement explicit,
we define a unitary operator

$\hat{\mathsf{P}}(t)=\hat{\Pi}_{g}\hat{D}(\alpha_{g}(t))+\hat{\Pi}_{e}\hat{D}(\alpha_{e}(t))+\hat{\Pi}_{f}\hat{D}(\alpha_{f}(t)),$
(155)

where $\hat{\Pi}_{a}=\ket{a}\\!\bra{a}$ are the projection operator onto the
energy eigenstate $\ket{a}$ of the qutrit and $\hat{D}(\alpha_{a}(t))$ are the
displacement operators (see the Eq.(39) for the definition). Intuitively,
$\hat{\mathsf{P}}$ entangles each projection $\hat{\Pi}_{a}$ of the qutrit
with a displacement operator of the resonator such that if the qutrit is in an
energy eigenstate, the resonator coherent state will be displaced to the
vacuum state. For the subsequent derivation, we follow the notation used in
[12] and use

$\hat{O}^{\mathsf{P}}=\mathsf{P}^{\dagger}\hat{O}^{\mathsf{P}}\mathsf{P}$
(156)

to denote any operator $\hat{O}$ in the displaced frame.

In the new frame, the density operator of the composite system is given by

$\hat{\rho}^{\mathsf{P}}(t)=\mathsf{P}^{\dagger}\hat{\rho}(t)\mathsf{P},$
(157)

where, to simplify the notation, to will start to use
$\hat{\rho}=\hat{\rho}_{\mathcal{SR}}$. In addition, if we define

$\hat{\rho}^{\mathsf{P}}_{nmab}(t)=\bra{n,a}\hat{\rho}^{\mathsf{P}}(t)\ket{m,b}$
(158)

to be the matrix element of $\hat{\rho}^{\mathsf{P}}$ in the energy basis of
the qutrit and the number basis of the resonator, then

$\hat{\rho}^{\mathsf{P}}=\sum_{n,m=0}^{\infty}\sum_{a,b\in\\{g,e,f\\}}\hat{\rho}^{\mathsf{P}}_{nmab}\ket{n,a}\\!\bra{m,b}.$
(159)

Our goal is to find the time evolution of the qutrit reduced density operator,
i.e.,

$\hat{\rho}_{\mathcal{S}}(t)=\Tr_{\mathcal{R}}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\rho}(t)\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}=\Tr_{\mathcal{R}}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\mathsf{P}\hat{\rho}^{\mathsf{P}}(t)\mathsf{P}^{\dagger}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}.$
(160)

By using Eq.(155), we obtain

$\displaystyle\hat{\rho}_{\mathcal{S}}(t)$
$\displaystyle=\sum_{n}\Big{(}\rho^{\mathsf{P}}_{nngg}\ket{g}\\!\bra{g}+\rho^{\mathsf{P}}_{nnee}\ket{e}\\!\bra{e}+\rho^{\mathsf{P}}_{nnff}\ket{f}\\!\bra{f}\Big{)}$
$\displaystyle\ \ \ \
+\sum_{n,m}\Big{(}\lambda_{nmmn}^{ge}\ket{g}\\!\bra{e}+\lambda_{mnnm}^{ge*}\ket{e}\\!\bra{g}\Big{)}$
$\displaystyle\ \ \ \
+\sum_{n,m}\Big{(}\lambda_{nmmn}^{gf}\ket{g}\\!\bra{f}+\lambda_{mnnm}^{gf*}\ket{f}\\!\bra{g}\Big{)}$
$\displaystyle\ \ \ \
+\sum_{n,m}\Big{(}\lambda_{nmmn}^{ef}\ket{e}\\!\bra{f}+\lambda_{mnnm}^{ef*}\ket{f}\\!\bra{e}\Big{)},$
(161)

where

$\displaystyle\lambda_{nmpq}^{ge}(t)$
$\displaystyle=\rho^{\mathsf{P}}_{nmge}e^{-\mathrm{i}\imaginary(\alpha_{e}\alpha_{g}^{*})}d_{pq},$
(162) $\displaystyle\lambda_{nmpq}^{gf}(t)$
$\displaystyle=\rho^{\mathsf{P}}_{nmgf}e^{-\mathrm{i}\imaginary(\alpha_{f}\alpha_{g}^{*})}d_{pq},$
(163) $\displaystyle\lambda_{nmpq}^{ef}(t)$
$\displaystyle=\rho^{\mathsf{P}}_{nmef}e^{-\mathrm{i}\imaginary(\alpha_{f}\alpha_{e}^{*})}d_{pq},$
(164)

with $d_{pq}(t)=\bra{p}\hat{D}(\beta_{ge})\ket{q}$. To arrive at Eq.(161), we
have used the fact that

$\sum_{p}\bra{m}\hat{D}^{\dagger}(\alpha)\ket{p}\bra{p}\hat{D}(\alpha)\ket{n}=\delta_{mn}$
(165)

since $\hat{D}(\alpha)$ is unitary. Once the matrix elements of
$\rho^{\mathsf{P}}$ (and thus $\lambda_{nmpq}^{ab}$) are known, the matrix
elements of the qutrit reduced density operator can be computed trivially. For
example,

$\hat{\rho}_{\mathcal{S},gg}=\bra{g}\hat{\rho}_{\mathcal{S}}\ket{g}=\sum_{n}\rho^{\mathsf{P}}_{nngg}$
(166)

and

$\hat{\rho}_{\mathcal{S},ge}=\bra{g}\hat{\rho}_{\mathcal{S}}\ket{e}=\sum_{n,m}\lambda_{nmmn}^{ge}.$
(167)

###  B.2 Master Equation in the Displaced Frame

The density operator in the displaced frame satisfies the master equation

$\displaystyle\dot{\hat{\rho}}^{\mathsf{P}}$
$\displaystyle=-\frac{\mathrm{i}}{\hbar}\Big{[}\hat{H}_{\text{eff}}^{\mathsf{P}},\hat{\rho}^{\mathsf{P}}\Big{]}-\hat{\mathsf{P}}^{\dagger}\dot{\hat{\mathsf{P}}}\hat{\rho}^{\mathsf{P}}-\hat{\rho}^{\mathsf{P}}\dot{\hat{\mathsf{P}}}^{\dagger}\hat{\mathsf{P}}+\kappa\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{a}^{\mathsf{P}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}^{\mathsf{P}}$
$\displaystyle\ \ \ \
+\gamma_{1,ge}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ge}^{\mathsf{P}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}^{\mathsf{P}}+\gamma_{1,gf}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{gf}^{\mathsf{P}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}^{\mathsf{P}}$
$\displaystyle\ \ \ \
+\gamma_{1,ef}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{ef}^{\mathsf{P}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}^{\mathsf{P}}+\frac{\gamma_{\phi,ge}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ge}^{\mathsf{P}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}^{\mathsf{P}}$
$\displaystyle\ \ \ \
+\frac{\gamma_{\phi,gf}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,gf}^{\mathsf{P}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}^{\mathsf{P}}+\frac{\gamma_{\phi,ef}}{2}\mathcalboondox{D}\mathopen{\scalerel*[5.5pt]{\Big{[}}{\ensurestackMath{\addstackgap[1.5pt]{\big{[}}}}}\hat{\sigma}_{z,ef}^{\mathsf{P}}\mathclose{\scalerel*[5.5pt]{\Big{]}}{\ensurestackMath{\addstackgap[1.5pt]{\big{]}}}}}\hat{\rho}^{\mathsf{P}}.$
(168)

As for any time-dependent unitary transformation, the extra terms
$-\hat{\mathsf{P}}^{\dagger}\dot{\hat{\mathsf{P}}}\hat{\rho}^{\mathsf{P}}-\hat{\rho}^{\mathsf{P}}\dot{\hat{\mathsf{P}}}^{\dagger}\hat{\mathsf{P}}$
appear in the new master equation to eliminate the readout drive terms in the
original master equation. Moreover, the Hamiltonian of the combined system
(qutrit + cavity) still takes the form

$\displaystyle\hat{H}_{\text{eff}}/\hbar$
$\displaystyle=\hat{H}_{\mathcal{SR},\text{rot}}^{\text{disp}}/\hbar$
$\displaystyle=\tilde{\omega}_{\text{q}}\hat{\Pi}_{e}+(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\hat{\Pi}_{f}+\Delta_{\text{rd}}\hat{a}^{\dagger}\hat{a}$
$\displaystyle\ \ \ \
+\chi_{\text{qr}}\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{\Pi}_{e}+2\hat{\Pi}_{f}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}\hat{a}^{\dagger}\hat{a}-\left(\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a}\right).$
(169)

Note, however, we have kept the qutrit decay (i.e., $\gamma_{1,ab}$) in the
master equation for full generality.

To begin simplifying each term, we will need various operators rewritten in
the displaced frame. For the cavity operators, we have

$\hat{a}^{\mathsf{P}}=\hat{a}+\left(\alpha_{g}\hat{\Pi}_{g}+\alpha_{g}\hat{\Pi}_{g}+\alpha_{g}\hat{\Pi}_{g}\right)=\hat{a}+\hat{\Pi}_{\alpha},$
(170)
$\displaystyle\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{a}^{\dagger}\hat{a}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}^{\mathsf{P}}$
$\displaystyle=\hat{a}^{\dagger}\hat{a}+\hat{a}^{\dagger}\hat{\Pi}_{\alpha}+\hat{a}\hat{\Pi}_{\alpha}^{\dagger}+\hat{\Pi}_{\alpha}^{\dagger}\hat{\Pi}_{\alpha}$
$\displaystyle=\hat{a}^{\dagger}\hat{a}+\hat{a}^{\dagger}\hat{\Pi}_{\alpha}+\hat{a}\hat{\Pi}_{\alpha}^{\dagger}$
$\displaystyle\ \ \ \
+|\alpha_{g}|^{2}\hat{\Pi}_{g}+|\alpha_{e}|^{2}\hat{\Pi}_{e}+|\alpha_{f}|^{2}\hat{\Pi}_{f},$
(171)

where we denote

$\hat{\Pi}_{\alpha}(t)=\alpha_{g}(t)\hat{\Pi}_{g}+\alpha_{g}(t)\hat{\Pi}_{g}+\alpha_{g}(t)\hat{\Pi}_{g}.$
(172)

Similarly, the operators associated with the qutrit subspace in the displaced
frame are given by

$\hat{\sigma}_{z,ab}^{\mathsf{P}}=\hat{\sigma}_{z,ab}\ \ \text{ and }\ \
\hat{\sigma}_{ab}^{\mathsf{P}}=\hat{\sigma}_{ab}\hat{D}^{\dagger}(\alpha_{a})\hat{D}(\alpha_{b}).$
(173)

First, we start with the transformed Hamiltonian. By using Eq.(170)-(173), we
obtain

$\displaystyle\frac{\hat{H}_{\text{eff}}^{\mathsf{P}}}{\hbar}$
$\displaystyle=\tilde{\omega}_{\text{q}}\ket{e}\\!\bra{e}+(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\ket{f}\\!\bra{f}$
$\displaystyle\ \ \ \
-\left[\epsilon\Big{(}\hat{a}^{\dagger}+\hat{\Pi}_{\alpha}^{\dagger}\Big{)}+\epsilon^{*}\Big{(}\hat{a}+\hat{\Pi}_{\alpha}\Big{)}\right]$
$\displaystyle\ \ \ \
+\Big{[}\Delta_{\text{rd}}+\chi_{\text{qr}}\mathopen{\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}}\hat{\Pi}_{e}+2\hat{\Pi}_{f}\mathclose{\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}}\Big{]}$
$\displaystyle\ \ \ \ \ \ \ \
\times\Big{(}\hat{a}^{\dagger}\hat{a}+\hat{a}^{\dagger}\hat{\Pi}_{\alpha}+\hat{a}\hat{\Pi}_{\alpha}^{\dagger}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
+|\alpha_{g}|^{2}\hat{\Pi}_{g}+|\alpha_{e}|^{2}\hat{\Pi}_{e}+|\alpha_{f}|^{2}\hat{\Pi}_{f}\Big{)}$
$\displaystyle=\tilde{\omega}_{\text{q}}\ket{e}\\!\bra{e}+(2\tilde{\omega}_{\text{q}}+\alpha_{\text{q}})\ket{f}\\!\bra{f}$
$\displaystyle\ \ \ \
-\Big{(}\epsilon\hat{a}^{\dagger}+\epsilon^{*}\hat{a}\Big{)}-\Big{(}\epsilon\hat{\Pi}_{\alpha}^{\dagger}+\epsilon^{*}\hat{\Pi}_{\alpha}\Big{)}$
$\displaystyle\ \ \ \